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ABSTR?^CT 

A  general  combustion  and  heat  transfer  model  for  porous 
media  subject  to  Darcy  flow  is  formulated.   The  transient 
one  dimensional  model  treats  the  combustion  process  in  two 
phases.   During  the  initial  phase,  combustion  occurs  within 
the  porous  medium.   The  second  phase  occurs  when  the  exo- 
thermic reaction  moves  to  the  air  inlet  surface  of  the  medium 
resulting  in  surface  recession.   The  temperature  dependency 
of  the  system  parameters  and  thermophysical  properties  is 
taken  into  account.   An  analysis  of  combustion  in  a  carbon 
porous  medium  is  presented,  as  well  as  an  assessment  of  the 
accuracy  of  the  model. 
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I.   INTRODUCTION 

Studies  were  conducted  at  the  Naval  Weapons  Center,  China 
Lake,  California  to  assess  the  survivability  of  the  F-18 
aircraft  when  exposed  to  open  pool  fires.   A  portion  of  the 
F-18  aircraft  is  constructed  of  graphite-epoxy  materials. 
When  this  composite  material  is  exposed  to  a  severe  thermal 
environment,  the  epoxy  burns  off  at  a  low  temperature  (about 
400-500  degrees  Fahrenheit)  leaving  behind  a  porous  graphite 
mat.   Experiments  were  conducted  by  J.  Fontenot  [1]  to  evaluate 
the  combustion  behavior  of  the  porous  graphite  or  carbon  mat. 
Due  to  the  complex  behavior  of  the  combustion  observed  for  the 
carbon  fibers,  a  mathematical  model  to  simulate  this  behavior 
was  developed  by  Vatikiotis  [2] .   The  model  was  formulated  to 
consider  a  porous  medium  composed  of  cylindrically  shaped 
fibers,  typical  of  fiber  reinforced  composite  materials.   The 
details  of  the  initial  effort  were  presented  by  Vatikiotis. 
Significant  improvements  have  been  made  to  the  initial  combus- 
tion model,  and  the  model  has  been  extended  to  consider 
spherically  shaped  particles,  typical  of  thermal  flow  reactors. 
The  initial  work  of  Vatikiotis  [2],  where  it  applies  to  the 
present  investigation  is  included  in  the  discussion  for  the 
sake  of  completeness.   Other  analytical  investigations  taking 
a  similar  approach,  but  concerned  with  other  aspects  of  com- 
bustion or  thermally  active  porous  media,  have  been  performed 
by  Kordylewski  [3] ,  Sahota  and  Pagni  [4] ,  and  Mehta,  Sams  and 
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Luss  [5] .   Kordylewski  performed  a  numerical  analysis  of  homo- 
geneous combustion  in  a  two  dimensional  reactor.   The  objec- 
tive of  Kordylewski ' s  investigation  was  to  show  that  the  flow 
within  the  reactor  influenced  the  conditions  necessary  for 
ignition.   Sahota  and  Pagni  considered  the  transient  behavior 
of  heat  transfer  in  a  porous  medium  when  exposed  to  a  fire. 
Their  analysis  was  part  of  an  overall  investigation  for  deter- 
mining the  internal  stresses  in  a  structure  caused  by  a  fire. 
Combustion  of  the  porous  medium  was  not  a  consideration. 
Mehta,  Sams,  and  Luss  investigated  the  transient  temperature 
response  of  a  packed-bed  reactor  when  the  entering  fluid 
temperature  was  decreased.   Their  model  assumed  that  the 
temperatures  of  the  porous  solid  and  the  fluid  were  the  same. 

Slattery  [6],  Yaron  [7],  and  Kassoy  [8]  present  three 
approaches  for  mathematically  modelling  a  porous  media. 
Slattery' s  and  Yaron ' s  approaches  were  similar  in  using  inte- 
gral methods  in  their  development.   Slattery  considered  only 
mass  transfer  whereas  Yaron  treated  transport  processes  in 
general.   Kassoy 's  approach  was  based  on  performing  Taylor 
series  expansions  of  energy  transport  processes  with  respect 
to  a  differential  volume  of  a  porous  medium.   An  energy  balance 
yielded  two  differential  heat  transfer  equations,  one  for  the 
porous  solid  and  one  for  the  fluid.   A  fundamental  requirement 
of  this  approach  is  that  the  particle  diameters  remain  small 
with  respect  to  the  thickness  of  the  porous  medium. 

Although  the  above  investigations  were  concerned  with 
specific  aspects  of  thermally  active  porous  media,  the  present 
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combustion  model  was  developed  to  treat  a  larger  class  of 
problems  (i.e.,  combined  mass  transfer,  heat  transfer,  and 
combustion) .   In  contrast  to  other  investigations,  the  present 
model  takes  into  account  the  temperature  dependency  of  all 
the  thermophysical  properties  of  the  air  and  the  geometric 
parameters  of  the  porous  medium.   The  physical  parameters  such 
as  porosity,  permeability,  and  thickness  affect  the  magnitude 
and  behavior  of  pressure-driven  air  flow  within  the  porous 
medium.   As  discussed  by  Kordylewski  [3] ,  the  air  flow  influ- 
ences the  conditions  for  ignition.   Therefore,  sensitivity 
analyses  were  performed  using  the  combustion  model  showing 
the  effects  of  these  parameters  on  combustion.   In  addition, 
a  functional  relation  between  the  internal  pore  velocity  and 
an  initial  uniform  temperature  distribution  needed  to  sustain 
combustion  is  presented.   The  initial  work  of  Vatikiotis  [2] 
showed  that  the  initial  conditions  are  important  in  deter- 
mining whether  or  not  sustained  combustion  will  occur.   An 
analysis  is  presented  showing  the  dependence  of  combustion 
on  the  shape  of  the  initial  temperature.   In  addition,  the 
effects  of  the  boundary  conditions  are  examined  using  the 
combustion  model.   Specifically,  there  is  a  change  in  the 
behavior  of  combustion  when  insulated  boundaries  on  the  porous 
solid  are  changed  to  permit  heat  transfer  by  radiation.   The 
Semenov  model  is  used  to  explain  the  behavior  observed  during 
the  analyses  discussed  above.   Although  the  Semenov  model 
was  originally  proposed  for  homogeneous  combustion,  the 
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results  show  that  the  Semenov  model  can  be  applied  to  hetero- 
geneous combustion  in  porous  media. 
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II.   DESCRIPTION  OF  THE  PROBLEM 

The  problem  under  investigation  is  that  of  a  porous  slab, 
with  a  pressure  differential  across  its  thickness,  in  a  severe 
thermal  environment.   The  effect  of  the  pressure  differential 
is  to  induce  air  flow  through  the  porous  medium.   This  internal 
air  flow  produces  two  opposing  effects,  (1)  internal  convection 
heat  transfer,  and  (2)  a  supply  of  oxygen  for  heat  generation 
by  combustion.   Whether  the  porous  medium  moves  towards  sus- 
tained combustion  or  extinguishment  depends  on  the  interaction 
of  these  effects.   A  mathematical  model  was  formulated  to 
provide  an  understanding  of  this  interaction  and  its  effect 
on  thermal  behavior. 

The  mathematical  model  was  formulated  as  follows.   Energy 
balances  on  the  porous  solid  and  the  convected  air  provide 
heat  transfer  equations  for  each.   The  heat  transfer  mechan- 
isms included  in  the  model  are  (1)  conduction,  (2)  convection, 
and  (3)  radiation.   In  addition,  nonvolatile  combustion  is 
included  in  the  porous  solid  heat  transfer  equation  as  a  heat 
generation  term  of  Arrhenius  type. 

The  conservation  of  species  law  was  applied  to  the  oxygen 
molecule  concentration.   The  resulting  mass  transfer  equation 
includes  the  transfer  mechanisms  of  (1)  molecular  diffusion, 
and  (2)  species  transport  by  convection.   A  teim  accounting 
for  the  consiomption  of  oxygen  due  to  combustion  is  included 
as  well.   Darcy's  law  was  taken  as  the  constitutive  equation 
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for  flow  through  the  porous  medium.  The  pressure  gradient, 
needed  to  calculate  pore  velocity,  is  obtained  by  solving  a 
combined  Darcy's  law  and  continuity  equation. 

As  the  behavior  of  the  system  occurs  over  a  large  tempera- 
ture range,  all  thermophysical  properties  which  depend  upon 
temperature,  such  as  conductivity,  viscosity,  and  air  density, 
are  treated  as  temperature  dependent  properties.   Moreover, 
as  combustion  consumes  the  porous  solid,  changes  in  porosity 
and  permeability  are  accounted  for  in  the  transient  analysis. 

In  order  to  provide  a  general  model,  a  number  of  boundary 
conditions  are  included  in  the  formulation  of  the  model .   A 
detailed  discussion  of  boundary  conditions  is  presented  in 
Section  III.H.   The  final  system  of  equations  to  be  solved  are 
four  transient,  nonlinear,  coupled  partial  differential  equa- 
tions.  The  three  heat  and  mass  transfer  equations  were  solved 
by  a  Galerkin  formulation  of  the  finite  element  method.   The 
remaining  combined  Darcy's  law-continuity  equation  was  solved 
by  a  shooting  method.   The  details  of  the  solution  procedure 
are  presented  in  Appendix  C. 
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III.   THEORY  AND  BACKGROUND 

A.   DESCRIPTION  OF  THE  POROUS  MEDIUM 

For  this  investigation,  a  porous  medium  is  defined  as  a 
solid  containing  interconnected  pores  (i.e.,  path  for  airflow). 
There  are  two  classes  of  porous  media,  consolidated  and  uncon- 
solidated.  Consolidated  porous  media  are  those  where  the 
solid  constituent  or  matrix  remains  rigid.   Consolidated  por- 
ous media  may  be  comprised  of  either  individual  particles 
attached  to  one  another  (e.g.,  sintered  metal),  or  a  solid 
where  a  portion  has  been  removed  (e.g.,  charcoal).   Unconsoli- 
dated porous  media  are  comprised  of  discrete  particles  (e.g., 
granular  beds) .   The  properties  of  an  unconsolidated  porous 
medium  can  change  if  the  porous  medium  is  agitated  or  settling 
takes  place.   Both  consolidated  and  unconsolidated  porous 
media  can  have  spatially  varying  properties . 

The  physical  properties  which  characterize  porous  media 
(i.e.,  allow  comparison  of  porous  media  without  regard  to  class 
or  form)  are  (1)  porosity,  (2)  specific  internal  area,  (3)  pore 
size  or  diameter,  and  (4)  tortuosity.   These  properties  are 
common  to  both  consolidated  and  unconsolidated  porous  media. 
Porosity,  p,  is  defined  as  the  ratio  of  void  volume  to  total 
volume.   The  specific  internal  area,  z,  is  the  ratio  of  in- 
ternal surface  area  to  bulk  volume.   The  dimension  is  the 
reciprocal  of  length.   Pore  diameter  is  used  to  describe  the 
pore  system  of  a  porous  medium.   The  actual  shape  and  size  of 
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a  pore  may  be  irregular  and  complex,  and  would  be  difficult 
to  describe  geometrically.   In  practice,  there  are  many  con- 
ventions used  to  define  the  pore  diameter.   For  this  inves- 
tigation, a  hydraulic  diameter  theory  is  used  and  is  discussed 
later  in  this  section.   The  tortuosity  of  a  porous  medium  is 
defined  as  the  ratio  of  the  length  of  the  flow  path  of  a  fluid 
particle  to  the  straight  line  distance.   Originally,  tortu- 
osity was  considered  a  kinematic  property,  but  was  subsequently 
adopted  for  characterizing  porous  media.   Tortuosity  is  also 
discussed  later  in  this  section.   Scheidegger  [9]  discusses 
various  methods  used  to  measure  properties  of  porous  media. 
Though  most  methods  discussed  are  based  on  experimental  tech- 
niques, the  above  properties  for  this  investigation  are  calcu- 
lated based  on  an  idealization  of  the  geometry  of  a  porous 
medium. 

The  porous  medium  was  modelled  as  shown  in  Figure  III.l. 
In  Figure  III.l,  D  is  the  distance  between  the  fiber  or  parti- 
cle centers,  and  d  is  the  particle  diameter.   The  geometry  is 
idealized  since  the  actual  porous  medium  may  be  irregular  in 
particle  diameter  and  distribution.   For  cylindrical  fibers, 
the  porosity,  which  was  defined  as  the  ratio  of  void  per  unit 
volume  is. 


=   1  -  j(d/D)^  (III.l) 


and  for  spherical  particles. 
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p   =   1  -  J(d/D)^  (III. 2) 


The  pore  diameter  is  obtained  by  an  expression  proposed  by 
Carman  [10] , 

6   =   4p/z  (III. 3) 

where  z  is  the  specific  internal  area  or  the  particle  surface 
area  per  unit  volume  of  porous  medium.   Expression  III. 3  is 
analogous  to  the  more  familiar  form  of  mean  hydraulic  diam- 
eter, 4V/A,  where  V  is  the  void  volume  and  A  is  the  wetted 
surface  area.   The  specific  internal  area,  z,  based  on  the 
idealized  geometry  of  Figure  III.l  is  calculated  by. 


z   =   i-TTd/D^  (III. 4) 


for  cylindrical  fibers,  and  by. 


z   =  ^ird^/D^  (III. 5: 


for  spherical  particles.   Expressions  III.  4  and  III. 5  assume 
one  half  the  total  internal  surface  area  is  effective  for 
convection  heat  transfer.   This  fractional  amount  of  total 
area  was  an  estimate  based  on  Fontenot's  experimental  results 
and  does  not  apply  to  porous  media  in  general.   A  more  general 
approach,  as  presented  by  Scheidegger  [9] ,  is  the  Kozeny 
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equation  given  by. 


z   =  isp^/m)-^^'^  (III. 6) 


In  equation  III.  6,  m  is  the  specific  permeability  or  hydraulic 
conductivity  (discussed  in  Section  III.B),  and  s  is  the  Kozeny 
constant.   The  Kozeny  constant  is  usually  taken  as  .2  for 
specific  surface  calculations.   Advantages  to  using  the  Kozeny 
equation  are  (1)  the  calculated  values  of  specific  internal 
area  are  in  fair  agreement  with  the  experimental  values,  and 
(2)  the  calculation  is  independent  of  particle  shape.   A 
disadvantage  is  that  the  Kozeny  equation  fails   for  highly 
porous  fibrous  media.   The  tortuosity,  t,  is  the  ratio  of  the 
length  of  the  flow  path  for  a  fluid  particle  to  the  straight 
line  distance.   For  the  geometric  configuration  of  Figure  III.l, 
the  tortuosity  depends  on  the  ratio  d/D.   Carman  [10]  pre- 
sents a  table  of  measured  tortuosity  factors  of  various  materials 
and  geometries ,  and  points  out  the  differences  between  the 
analytical  determinations  of  tortuosity.   Here,  we  adopt  his 
recommended  value  of  1.4.   Since  particle  size  decreases  as 
the  carbon  is  consumed,  all  geometric  properties  which  depend 
on  fiber  or  particle  diameter  are  functions  of  time  and  posi- 
tion.  The  model  assumes  that  the  carbon  matrix  remains  rigid 
as  the  particle  diameter  decreases,  and  thus,  porosity  in- 
creases with  combustion.   In  his  discussion  on  the  packing 
theory  of  spheres,  Scheidegger  [9]  states  the  highest  reported 
porosity  for  spherical  particles  in  a  stable  configuration  is 
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.875.   Carman  [10]  reports  of  investigations  of  fibrous  porous 
media  with  porosities  as  high  as  .99.   The  change  in  particle 
diameter  is  accounted  for,  and  will  be  discussed  in  the  sec- 
tion on  the  Arrhenius  expression  (Section  III.D)  for  carbon 
combustion. 

B.   DARCY'S  LAW  AND  PORE  VELOCITY 

Reynolds  number  for  porous  media  is  defined  by 


Re  =   (p  ud)/y  (III. 7) 

a 


where  u  is  the  local  pore  velocity,  p   is  the  mass  density  of 
air,  and  ]i   is  the  dynamic  viscosity.   Depending  on  the  magni- 
tude of  the  Reynolds  number,  the  motion  of  the  fluid  may  be 
dominated  by  molecular,  viscous,  or  inertial  effects.   Most 
investigations  of  fluid  flow  in  porous  media  have  been  in  the 
range  of  Reynolds  number  where  the  flow  is  dominated  by  vis- 
cous and  inertial  effects.   The  Navier-Stokes  equations  are 
applicable  in  describing  the  fluid  motion  in  this  range  of 
Reynolds  number.   However,  because  of  the  convoluted  geometry 
and  the  necessary  boundary  conditions  (i.e.,  u  =  0  at  a  solid- 
fluid  interface) ,  solution  of  the  Navier-Stokes  equations  are 
difficult  for  a  porous  medium.   Extensive  experimental  work, 
as  discussed  by  Scheidegger  [9]  has  shown  that  fluid  flow  in 
porous  media  is  governed  by  Darcy ' s  law  for  the  range  of 
Reynolds  number  where  viscous  effects  dominate.   The  upper 
limit  (high  velocity  end)  Reynolds  number  reported  in  the 
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experimental  investigations  varied  from  .1  to  75.   As  Reynolds 
number  increases,  inertial  effects  increase.   However,  Boffa 
[11]  has  shown  that  for  a  fixed  Reynolds  number,  inertial 
effects  diminish  with  increasing  air  temperature .   For  the 
problems  covered  in  this  investigation,  the  Reynolds  number 
did  not  exceed  a  value  of  5. 

Neglecting  body  forces,  Darcy's  law  for  one  dimensional 
flow  is. 


Q   =   -  -(^)  (III. 8) 

^       y  dx 


where  Q  is  the  filter  velocity  or  the  volumetric  flow  rate  per 
unit  of  cross-sectional  area,  and  dP/dx  is  the  pressure  gradi- 
ent.  The  specific  permeability  or  the  hydraulic  conductivity 
of  the  porous  medium  is  defined  by. 


m  =   96  p  (6/t)^  cm. 9) 


Expression  III. 9  is  based  on  a  capillaric-serial  model  given 
by  Scheidegger  [9] .   It  should  be  noted  that,  in  a  physical 
sense,  permeability  and  porosity  are  not  related.   Porosity 
is  a  quantifiable  property  of  the  porous  medium.   Whereas, 
permeability  is  a  constitutive  property  (i.e.,  a  property 
specified  by  a  constitutive  equation,  Darcy's  law).   As  in 
the  case  of  Fourier's  law  of  heat  transfer,  Darcy's  law  is 
not  derived  from  first  principles,  but  has  been  obtained 
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through  exhaustive  experimental  analyses.   The  Dupuit- 
Forcheimer  assumption,  presented  by  Carman  [10] ,  relates 
the  local  pore  velocity  to  the  filter  velocity  by, 

Q   =   pu  (III. 10) 

The  hypothesis  of  the  Dupuit-Forcheimer  assumption  is  that  the 
local  pore  velocity  is  greater  than  the  filter  velocity. 
Noting  that  the  actual  velocity  in  a  single  pore  is  not  a 
constant,  but  rather  a  function  of  the  location  within  the 
pore,  the  Dupuit-Forcheimer  assumption  defines  an  "average" 
velocity  within  the  pore. 

The  continuity  equation  for  one  dimensional  fluid  flow  in 
porous  media  with  nonconstant  porosity  distribution  is. 


^t^  +  PP^  1$  =   0  (III. 11) 


Substituting  in  the  Dupuit-Forcheimer  assumption  and  Darcy ' s 
Law,  the  continuity  equation  becomes. 


^  ^  ^p^  3x  ^  m  33^   U  d^^dSl       ^   dt  -      ^    (III. 12) 


Equation  III. 12  along  with  the  boundary  conditions  can  be 
integrated  to  provide  the  pressure  and  the  pressure  gradient 
distributions  through  the  porous  medium.   The  boundary  condi- 
tions are  P(0)  =  P  ,  the  ambient  pressure,  and  P (L)  =  P^ .   The 

00  ■^  Xj 
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pressure  differential  across  the  porous  medium  is  AP  =  P   -  P^ 
Knowing  the  pressure  gradient  distribution,  the  pore  velocity 
distribution  is  obtained  from  Darcy's  law  and  the  Dupuit- 
Forcheimer  assumption.   Derivations  of  equation  III. 12  and 
the  continuity  equation  are  presented  in  Appendix  A  and  B, 
respectively . 

C.   SEMENOV  MODEL  OF  COMBUSTION 

In  this  work  we  adopt  the  combustion  model  of  N.   N.  Semenov 
as  described  in  the  texts  of  Frank- Kamenetskii  [12]  and  Vulis 
[13] .   A  brief  discussion  of  those  features  of  the  model  which 
relate  to  the  present  investigation  will  be  given.   Fundamental 
to  the  model  is  the  relation  of  reaction  rate  to  temperature, 
and  the  interaction  of  heat  generation  and  heat  transfer. 
The  reaction  rate,  R  ,  is  represented  by  the  Arrhenius  law  for 
a  simple  n-th  order  reaction. 


R   =   A^^expC-^)  (III.  13) 

^  R  T 

u 

where  A   is  the  characteristic  time  of  the  chemical  reaction, 
E  is  the  activation  energy,  R   is  the  universal  gas  constant, 
T  is  the  absolute  temperature,  and  0  is  the  oxygen  concentra- 
tion.  A  simple  reaction  is  one  in  which  the  rate  depends  on 
the  concentrations  of  the  reactants  and  not  on  the  products. 
The  heat  generated  by  the  exothermal  reaction  is  obtained  by 
multiplying  R  by  the  heat  of  combustion. 
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In  Figure  III. 2,  heat  generation  R   is  plotted  versus 
temperature.   The  resulting  curve  is  referred  to  as  the  S- 
curve  due  to  its  sinuous  shape.   The  S-curves  are  distinguished 
by  two  regions.   In  region  I,  the  reaction  rate  and  tempera- 
ture are  low  and  as  a  result,  there  is  an  excess  supply  of 
oxygen.   In  this  region  the  reaction  is  controlled  by  the 
temperature,  and  is  called  the  kinetic  regime  of  combustion. 
In  the  kinetic  regime,  the  reaction  rate  increases  exponen- 
tially with  increasing  temperature.   In  region  II,  the  higher 
reaction  rates  are  only  slightly  affected  by  the  higher  tem- 
peratures, and  the  reaction  is  limited  by  the  availability  of 
oxygen.   Region  II  is  called  the  diffusion  regime  of  combus- 
tion.  It  should  be  noted  that  S-curves  are  not  obtained  by 
plotting  R  versus  temperature  for  a  constant  oxygen  concen- 
tration.  The  flattening  of  the  S-curve  at  the  higher  tempera- 
tures is  due  to  decreasing  oxygen  concentration  at  those 
temperatures.   The  S-curve,  which  represents  the  behavior  of 
an  actual  system,  demonstrates  the  strong  coupling  between 
temperature  and  oxygen  concentration. 

In  addition  to  the  S-curves,  Figure  III. 2  shows  two  con- 
vection heat  transfer  curves,  q^,  and  qn2'  ^'^  ^^^   flow  tem- 
peratures, T,  and  T^,    respectively.   The  equation  for  the 
heat  transfer  lines  is  of  the  form. 


q^   =   h(Tg  -  T^)  (III. 14) 
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where  h  is  a  convection  heat  transfer  coefficient,  and  T  and 

T  are  the  solid  surface  and  air  temperatures,  respectively, 
a 

The  intersection  of  the  heat  generation  and  heat  transfer 

curves,  R  =   qn  ,    defines  the  stable  quasi-stationary  point  A. 
g     X, 

For  a  system  in  thermal  equilibrium,  the  solid  and  air  tempera- 
tures will  remain  at  T   and  T, ,  respectively.   Perturbing  the 
system  from  one  equilibrium  position  to  another  causes  the 
quasi-stationary  point  to  adjust  accordingly  (e.g.,  raising 
the  air  flow  temperature  from  T-,  to  T^  will  cause  the  quasi- 
stationary  point  to  move  from  T,  to  T^) .   At  point  I,  the  R 
and  q^,  curves  are  tangent,  and  an  infinitesimal  increase  in 
temperature  will  result  in  a  jump  in  reaction  temperature  from 
Tj  to  T„.   At  point  I,  which  is  defined  as  the  "critical  igni- 
tion condition",  the  reaction  moves  from  the  kinetic  regime 
to  the  diffusion  regime  of  combustion.   Temperature  T^  is 
referred  to  as  the  ignition  temperature.   Although  the  des- 
cription of  the  process  just  presented  is  for  a  surface  in 
which  only  heat  transfer  by  convection  occurs,  the  underlying 
ideas  carry  over  for  the  total  heat  transfer  at  a  point  in  the 
porous  medium.   The  heat  transfer  can  include  contributions 
from  conduction  and  radiation,  as  well  as  from  convection. 
In  his  original  work,  Semenov  [14]  developed  his  theory 
to  explain  reaction  behavior  within  homoegenous  gas  mixtures. 
Frank-Kamenetskii  [12]  later  adapted  this  theory  to  describe 
the  behavior  associated  with  heterogeneous  mixtures  (i.e.,  a 
solid  surface  and  air) .   In  their  study  of  coal  combustion. 
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Thomas,  Stevenson,  and  Evans  [15]  state  that  a  "temperature 
jump",  as  described  by  Frank-Kamenetskii,  may  or  may  not 
occur  depending  on  the  partial  pressure  of  oxygen.   This  pres- 
ent investigation  does  not  consider  this  aspect  of  the  com- 
bustion theory.   However,  analyses  were  made  using  the 
present  model  to  determine  whether  the  reaction  takes  place 
in  the  kinetic  regime  or  in  the  diffusion  regime.   Moreover, 
the  Semenov  model  will  be  used  to  explain  these  results  in 
Chapter  V. 

D.   ARI^ENIUS  LAW  OF  REACTION  RATE 

For  carbon  reacting  in  air,  Parker  and  Hottel  [16]  pro- 
posed the  following  Arrhenius  expression  for  the  reaction 

2 
rate  in  units  of  kg-carbon/m  -s, 

^0 
R^   =   9.55  X  10^  —^   exp("'^'^^^^)  (III. 15) 

c         u  c 

where  P   is  the  partial  pressure  of  oxygen  in  atmosphere 
units,  and  R   is  the  universal  gas  constant  (1.986  cal/gmole- 
K) .   Equation  III. 15  assumes  a  simple  first  order  reaction 
for  the  combustion  of  carbon  in  air.   Frank-Kamenetskii  [12] 
has  shown  that  the  experimental  data  of  Parker  and  Hottel  is 
better  correlated  by  a  fractional  order  reaction,  where  the 
reaction  order,  n,  is  between  1/3  and  2/3.   Replotting  Parker's 
and  Hottel 's  data  for  n  =  1/2,  yields  the  reaction  rate  ex- 
pression used  in  the  present  analysis. 
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R   =   2.Q65xlQ^(R   (t))-*-^"^  expC^'^^^^)  (III. 16) 

^  "2  R   T 

u   c 

2 
In  expression  III. 16,  R  is  in  units  of  Ibm- carbon/ ft  - 

hr,  R_   is  the  gas  constant  for  oxygen  (48.29  ft-lbf/lbm-R) , 

2 

3 
(J)  is  in  Ibm/ft  ,  R   is  1.986  Btu/lbmole-R,  and  T   is  in  Rankines 

The  ideal  gas  law  was  used  to  replace  the  partial  pressure, 

P   .   Figure  III. 3  and  Table  III.l  show  the  results  of  plotting 

"2 
Parker's  and  Hottel ' s  data  for  values  of  n  equal  to  1,  2/3, 

1/2,  and  1/3.   The  values  presented  in  Table  III.l  are  in  S.I. 
units. 

In  order  to  determine  the  rate  of  heat  generation  and  the 
rate  of  oxygen  consumption,  the  chemical  reaction  for  the 
combustion  process  must  be  considered.   For  nonvolatile  com- 
bustion of  carbon  and  oxygen,  two  reactions  that  describe  the 
process  are. 


C  +  J  O   -^   CO  (III. 17) 


C  +  O2   -»■   CO2  (III. 18) 


The  ratio  of  the  mass  rates  of  carbon  monoxide  to  carbon 
dioxide  produced  increases  with  increasing  temperature. 
Arthur  [17]  presents  an  expression  for  the  rate  ratio  as  a 
function  of  temperature  (in  Kelvins) . 


CO/CO^   =   2500  expC^^^^)  (III. 19) 

^  T 

c 
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FIGURE  III. 3  Results  of  plotting  reaction  rate 

data  of  Parker  and  Hottel  [16]  for 
several  values  of  n. 
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R^  =  A  ^ , 

=      JV 


exp' 


^t 


gr-c 
cnrvs 


n 

A 

r  g-c-K"  1 

^         fgr-caf 

Cm-atm-sec 

L  mole 

1 
3 

1.19  X  10^ 

27500 

1 
2 

1.29  X  10^^ 

31800 

2 

3 

2.1   X  10^ 

37100 

1 

9.55X  10^ 

44000 

TABLE  III.l  Results  of  plotting  reaction  rate 

data  of  Parker  and  Hottel  [16]  for 
several  values  of  n. 
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As  a  result  of  this  dependence  on  temperature,  the  stoichio- 
metric ratio  and  the  heat  of  reaction  will  also  be  functions 
of  temperature.   Defining  the  fraction  of  carbon  monoxide 
being  produced  by. 


( CO/CO ^) 

^co  =  1  ^  (CO/CO,)  ^^^^-20) 


and  the  fraction  of  carbon  dioxide  as. 


^C02   "   1  +  (CO/CO2)  (III. 21) 


the  heat  of  combustion  is  then  expressed  as. 


^Hr  =   FcO  ^«C0  ^  ^C02  ^«C02  '"^•'2' 


Values  for  the  heats  of  combustion,  AH_-  and  AH^-  ,  as  func- 
tions of  temperature  were  obtained  from  the  JANAF  (Joint  Army, 
Navy,  and  Air  Force)  Tables  [18] .   The  stoichiometric  ratio 
(fuel  to  oxygen)  of  the  overall  reaction  is. 


^R  =  ^co  ^co/<fco2  ^co  ^  ^co  ^co^'       '"^•"' 


where  f^-  is  the  stoichiometric  ratio  for  reaction  III. 17  and 
fpQ   is  the  stoichiometric  ratio  for  reaction  III. 18.   The 
rate  of  heat  generation,  R  ,  and  the  rate  of  oxygen  consumption 


can  now  be  expressed  by. 
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R    =   AH„  R  (III. 24) 

g        Re 


\      =  ^'r     %  <"I-25) 


It  should  be  pointed  out  that  Parker's  and  Hottel's  work  [16] 
and  Arthur's  work  [17]  were  accomplished  using  specific  types 
of  carbon.   Smoot  and  Pratt  [19]  and  Frank-Kamenetskii  [12] 
list  tables  and  references  for  many  other  rate  expressions, 
depending  on  the  particular  type  of  carbon  used.   The  present 
model  was  developed  for  any  rate  expression  of  Arrhenius  type 
for  carbon  consumption.   For  this  investigation,  we  will  re- 
main with  Parker's  and  Hottel's  expression  III. 16  as  modified 
by  Frank-Kamenetskii . 

As  previously  stated,  the  particle  diameter  decreases  as 
combustion  progresses.   The  rate  of  decrease  depends  on  fhe 
amount  of  carbon  consumed  at  a  point  over  time.   Past  observa- 
tions have  shown  that  the  effect  is  significant  when  the  reac- 
tion is  concentrated  in  a  small  region  of  the  porous  medium. 
To  take  this  into  account,  expressions  for  the  time  rate  of 
change  of  the  diameter  as  a  function  of  reaction  rate  were 
derived.   The  equation  for  cylindrical  fibers  is. 


d   =   -2  R  zdV(ttp  d)  (III. 26) 


and  for  spherical  particles. 


d   =   -2  R  zD^/(7Tp  d^)  (III. 27) 
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where  p   is  the  bulk  mass  density  of  carbon.   The  diameter 
and  the  reaction  rate  are  functions  of  both  time  and  position. 


E.   HEAT  TRANSFER  EQUATIONS  FOR  POROUS  MEDIA 

In  a  previous  investigation  of  porous  media.  Green  and 
Perry  [20]  developed  heat  transfer  equations  for  the  solid 
and  fluid  phases.   The  model  included  the  three  basic  mechan- 
isms, (1)  physical  movement  of  the  fluid  with  its  own  heat 
content,  (2)  conduction  of  heat  through  the  solid  and  fluid 
phases,  and  (3)  convective  heat  transfer  between  the  solid 
and  the  fluid  phases.   Radiation  heat  transfer  was  assumed 
to  be  negligible,  and  combustion  was  not  a  consideration. 
The  differential  equations  for  the  solid  and  the  fluid  phases, 
obtained  from  energy  considerations, 

Ps=s<l-P'Tr  =  '^s<l-P>— T  *  '^^'^w  -  ^s>        (III. 28) 

3T  dT  3^T 

p  c  p-^     =   -pp  c  u-x-^  +  k  p ^   -  hz(T   -  T  )        (III.  29) 

^w  w^  3t       ^^w  w9x     w^2        w    s 

dX 

and  were  solved  numerically  by  Green  and  Perry  using  the 
finite  difference  method.   In  equations  III. 28  and  III. 29, 
c  is  specific  heat,  h  is  the  internal  convection  heat  trans- 
fer coefficient,  k  is  thermal  conductivity,  and  z  is  the  wetted 
surface  area  per  unit  volume.   The  s  and  w  subscripts  refer 
to  solid  and  fluid  properties,  respectively.   Their  model  does 
not  account  for  the  temperature  dependency  of  the  properties. 
Other  investigators  such  as  Sundaresan  [21] ,  Riaz  [22] , 
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Mendelsohn  [23] ,  and  Schneider  [24]  have  presented  a  single 
equation  approach  to  describe  heat  transfer  in  porous  media. 


[(l-p)c^p^+pc„p^l|f  =   -up„c„p||^  [(l-p)Vpk„jiq   (III. 30) 

oX 


Equation  III. 30  can  be  derived  from  equations  III. 28  and 
III. 29  by  assuming  the  solid  and  fluid  phase  temperatures  are 
equal .   Green  and  Perry  proposed  that  the  equal  temperature 
assumption  was  valid  for  systems  satisfying  the  following 
condition. 


hz  1/2  "^w 

^   =   (iT^)  -TT     >   -^42  (III. 31) 

K  p  U 
W*^ 


where  a   is  the  thermal  diffusivity  of  the  fluid. 
w 

In  the  present  investigation,  the  heat  transfer  equations 
are  generalized  to  also  include  (1)  radiation,  (2)  internal 
combustion,  (3)  temperature  dependency  of  properties,  and 
(4)  compressibility  effects  of  the  air. 

An  energy  balance  on  the  carbon  gives  the  heat  transfer 
equation. 


|3j((l-p)(k^+k^)i3£l-hz(T^-T^)+R  z   =   (l-p)p^c^^   (III. 32) 


where  k   is  a  pseudo  thermal  conductivity  to  account  for 
radiation  between  particles.   The  derivation  of  equation 
III. 32  is  presented  in  Appendix  A.   The  effective  conductivity, 
k  ,  of  the  porous  solid  was  proposed  by  Russel  [25] , 
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p.2/3,^„.p.2/3) 

k   =   k  r^ (III. 33) 

®      ^   ,2/3    ,  ^  a,,    ,2/3^  ,, 

where  k  and  k  are  the  bulk  thermal  conductivities  of  carbon 
c      a 

and  air,  respectively,  and  p'  =  (1-p) .   Russel's  expression, 
based  on  an  electrical  resistance  analogy,  can  be  used  for 
the  full  range  of  porosity,  from  0  to  1. 

The  pseudo  thermal  conductivity  due  to  radiation,  k  ,  is 
obtained  as  follows.   Treating  the  idealized  geometry  of  the 
porous  medium  as  a  series  of  closely  spaced  walls,  the  net 
radiation  heat  flux  between  two  of  the  walls  is. 


q^.   =   aeCrJ  -  T2)/(2-£)  (III. 34) 


where  £  is  the  emissivity,  and  a  is  the  Stef an-Boltzman  con- 
stant.   By  expanding  q   in  a  Taylor  series  about  T- /  an 
analogy  to  Fourier's  law  of  heat  transfer  by  conduction  yields 
the  effective  radiation  conductivity  as. 


k   =   4  a  £  5  tV(2-£)  (III. 35) 

r  c 


Expression  III. 35  is  then  used  for  k   in  equation  III. 32.   The 
details  of  the  above  derivation  are  presented  in  Appendix  B. 
The  application  of  expression  III. 35  to  porous  media  is  also 
proposed  by  Rohsenow  and  Hartnett  [26] .   In  a  more  rigorous 
development,  Whi taker  [27]  presents  an  alternative  approach 
for  treating  radiation  heat  transfer  in  porous  media. 
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Following  the  experimental  work  of  Yoshida,  Ramaswami,  and 
Hougen  [28] ,  the  internal  convection  heat  transfer  coefficient 
is  given  by  the  empirical  correlation. 


h   =   0.91  Re'~°*^^[i|^  c   G(c  u/k  )    ~'^^^]  (III. 36) 

a    a   a  rm 


where  ip    is  equal  to  .91  for  cylindrical  fibers,  and  equal  to 

1  for  spherical  particles,  G  is  a  pseudo  mass  velocity  given 

by  pp  u,  and  c   is  the  specific  heat  of  air  at  constant  pres- 
a        a 

sure.   The  fm  subscript  refers  to  the  properties  evaluated 
at  film  temperature.   Re'  is  a  pseudo  Reynolds  number  defined 
by. 

Re'   =   G/(zyi|;)  (III. 37) 

The  air  properties  vary  with  temperature,  and  h  varies  with 
temperature.   Finally,  the  reaction  rate  for  the  heat  genera- 
tion term  in  equation  III. 32  is  given  by  expression  III. 24. 

An  energy  balance  on  the  air  within  the  porous  medium 
provides  the  second  heat  transfer  equation  as, 

(III. 38) 

The  density  of  air,  p  ,  is  approximated  in  the  model  by  the 

a 

ideal  gas  law.   The  term,  3(pP)/3x,  is  due  to  the  compressi- 
bility of  the  air.   The  derivation  of  equation  III. 38  is 
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presented  in  Appendix  A.   All  properties  in  equation  III. 38 
vary  with  temperature.   The  properties  of  standard  air  were 
used  in  the  analysis.   During  combustion,  oxygen  in  the  air 
is  replaced  by  carbon  monoxide  and  carbon  dioxide.   Calcula- 
tions showed  that  the  gas  mixture  would  have  thermophysical 
properties  slightly  different  from  those  of  standard  air. 
In  a  "worst  case"  situation,  the  average  difference  between 
the  thermophysical  properties  of  air  and  a  mixture  of  .79  N2 
and  .21  CO2  is  7  percent.   This  mixture  assumes  that  the  oxy- 
gen is  totally  consumed.   Viscosity  had  the  largest  difference 
of  13  percent,  and  specific  heat  the  smallest  at  1.5  percent. 
The  properties  of  a  .79  N2  and  a  .21  CO  mijcture  are  approxi- 
mately those  of  air  (different  by  less  than  2  percent) .   At 
typically  observed  temperatures  (greater  than  1200  deg-F) ,  the 
combustion  gas  will  be  mostly  carbon  monoxide,  and  the  error 
introduced  by  assioming  the  properties  of  standard  air  is  mini- 
mum.  A  small  difference  in  the  properties  is  acceptable  since 
an  additional  molecule  mass  transfer  equation  for  either  carbon 
monoxide  or  carbon  dioxide  would  be  necessary  to  calculate 
the  effective  properties  of  the  gas  mixture.   The  methods  used 
to  obtain  the  thermophysical  properties  of  the  gas  mixtures 
above  were  those  of  Mason  and  Saxena  [29]  ,  and  Reynolds  and 
Perkins  [30]  .   The  polynomial  expressions  used  to  calculate 
the  thermophysical  properties  of  air  are  derived  in  Appendix  B. 

F.   OXYGEN  DIFFUSION  EQUATION  FOR  POROUS  MEDIA 

As  a  result  of  combustion,  the  heat  transfer  equations 
include  four  response  variables,  the  carbon  and  air  temperatures, 
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T  and  T  ,    the  concentration  of  oxygen,  4),  and  the  total 
c      a 

internal  pressure,  P.   The  fourth  field  equation  necessary 
to  complete  the  system  is  obtained  from  species  conservation 
considerations.   The  oxygen  molecule  transport  mechanisms 
included  in  the  model  are  (1)  molecular  diffusion  resulting 
from  concentration  gradients  (Pick's  Law),  (2)  convective 
mass  flow,  and  (3)  oxygen  consumption  due  to  combustion. 
Diffusion  resulting  from  pressure  and  temperature  gradients 
was  considered  negligible.   An  example  cited  by  Bird,  Stewart, 
and  Lightfoot  [21]  where  pressure  diffusion  is  important  is 
the  centrifuge  separator;  and  for  thermal  diffusion,  the 
Clusius-Dickel  column,  where  very  steep  temperature  gradients 
are  used  to  separate  complex  mixtures  of  organic  molecules. 
Thus,  the  oxygen  molecule  transfer  equation  for  a  porous 
medium  is. 


fe<P^el^)  -k<P"*'  -\^  =  Pff        ("1-35) 


The  derivation  of  equation  III. 39  is  presented  in  Appendix  A. 
The  second  order  term  in  equation  III. 39  results  from  Pick's 
law  of  diffusion.   The  effective  diffusivity  of  the  porous 
medium,  V    ,    is  a  function  of  pressure,  temperature,  and  pore 
geometry  of  the  medium.   Scheidegger  [9]  presents  several 
models  for  V      which  have  been  proposed  by  a  number  of  inves- 
tigators.  Briefly,  there  are  two  limiting  cases  of  diffusion, 
(1)  diffusion  associated  with  the  collision  between  gas  mole- 
cules, and  (2)  diffusion  associated  with  the  collision  of 
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gas  molecules  with  the  pore  walls  of  the  medium.   This  second 
phenomenon  is  discussed  by  Bennett  and  Myers  [32],  and  is 
referred  to  as  Knudsen  diffusion.   The  former  case  occurs 
when  the  mean  free  path  of  the  molecules  is  smaller  than  the 
pore  diameter.   The  expression  for  mean  free  path,  m ,    given 
by  Treybal  [33]  is. 


03   =   ^^(R^T/2tt  g^M)^/^  (III. 40) 


where  g   is  the  gravitational  constant  and  M  is  the  molecular 
weight.   From  equation  III. 40,  a  value  of  uj  =  5/100  is  obtained 
for  oxygen  for  the  geometry  of  typical  porous  media  encoun- 
tered in  the  analysis.   Thus,  collision  between  molecules  is 
the  governing  diffusion  mechanism.   A  semi-empirical  expres- 
sion, proposed  by  Gilliland  [34],  is  used  to  obtain  the  diffu- 
sion coefficient  of  oxygen  into  air.   Gilliland 's  expression 
is  given  by. 


V      =   435.7  T^/-^  (m"-^ +  m"-'-)/[P(V-^/^  + V^^)  ]         (III. 41) 
a     a    ^  o      ^      o 

2 
where  V    is  in  units  of  cm  /sec,  P  is  the  total  pressure  in 

Pa,  V  and  V-   are  the  molecular  volumes  of  air  and  oxygen, 
a      U2 

respectively.   M  and  M   are  the  molecular  weights  of  air  and 

a      U2 

oxygen.   The  values  of  V  and  V-   were  obtained  from  Holman 

a      O2 

[35]  as  29.9  and  7.4,  respectively.   The  effective  diffusivity 
proposed  by  Denbigh  and  Turner  [36]  for  a  porous  medium  is. 
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V        =  V/T  (III. 42) 

e 


Expression  III. 42  accounts  for  the  tortuous  path  the  oxygen 
molecules  follow  through  the  porous  medium.   Lastly,  the  oxy- 
gen consumption  term  in  equation  III.  39  is  given  by  expression 
III. 25. 

G.   THE  SURFACE  RECESSION  PROBLEM 

The  previous  system  of  equations  describes  the  combustion 
process  occurring  within  the  porous  medium.   As  combustion 
progresses,  the  oxygen  concentration  goes  to  zero,  and  as  a 
result,  the  reaction  moves  to  the  air  inlet  surface  of  the 
porous  medium.   Thereafter,  the  carbon  consumption  takes  place 
at  the  surface.   During  the  surface  recession  phase,  the  thick- 
ness of  the  porous  medium  decreases  with  time,  i.e.,  surface 
recession.   This  suggests  partitioning  the  problem  into  two 
parts,  (1)  initial  combustion  within  the  porous  medium,  and 
(2)  combustion  at  the  air  inlet  surface.   The  second  phase 
of  the  problem  is  formulated  as  a  moving  boundary  problem. 
The  governing  equations  for  the  moving  boundary  problem  are 
those  developed  previously  with  the  exception  of  the  oxygen 
molecule  mass  transfer  equation.   With  combustion  occurring 
at  the  surface,  oxygen  is  totally  consumed  in  a  shallow  region 
near  the  surface.   Experimentation  by  Koizumi  [37]  has  shown 
this  penetration  of  the  oxygen  to  be  1  to  2  particle  diameters. 
As  a  result,  the  oxygen  molecule  mass  transfer  equation  is 
not  needed  for  this  phase  of  the  problem,  and  the  heat  genera- 
tion can  be  treated  as  a  planar  source. 
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The  carbon  heat  transfer  equation  with  the  modified  heat 
generation  term  is. 


3T 

l"^|-[(1-P)  (k  +k  )^]-hz(T  -T  )+R'z6^(n  =  0) 
dT]  ^    e   r  3ri       c   a    g 


3T     3T 
(l-P^Pc^cfE8Tr^^Tt^l  (III. 43) 


where  ri  is  equal  to  x/L.   The  first  term  on  the  right  side  of 
equation  III. 43  arises  from  the  x  coordinate  becoming  a  func- 
tion of  time  during  surface  recession.   A  similar  term  appears 
in  the  air  temperature  and  combined  Darcy's  law-continuity 
equations  as  well.   Transformation  of  the  field  equations  from 
a  fixed  coordinate  to  a  moving  coordinate  system  is  presented 
in  Appendix  A.   As  was  observed  (see  Figure  V.2) ,  the  tempera- 
ture gradient  through  the  porous  medium  during  surface  reces- 
sion was  small.   As  a  result,  the  additional  term  may  be 
omitted  without  affecting  the  results.   The  moving  boundary 
problem  formulation  is  based  on  Crank's  [38]  extension  of  a 
method  proposed  by  Murray  and  Landis  [39]  .   The  thickness,  L, 
is  updated  continuously  by  evaluating  its  time  rate  of  change 
given  by. 


L  =   -(puf  )    (|)  /[(1-p)    p  ]  (III. 44 

x=0         x=0  ^ 


Expression  III. 44  is  obtained  from  a  mass  balance  on  the 
carbon  (i.e.,  the  time  rate  of  change  in  carbon  is  equal  to 
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the  carbon  consumption) .   Noting  that  all  the  oxygen  entering 
the  plate  is  consiimed  in  a  small  region  near  x  =  0,  the  carbon 
consumption  can  be  represented  in  terms  of  the  oxygen  flow 
and  the  stoichiometric  ratio  (i.e.,  R'  =  -puf_(J)  A  where  A 

• 

is  area  )  .   The  time  rate  of  change  for  carbon  is  (1-p) p  AL. 
The  planar  heat  source,  R' ,  in  equation  III. 43,  is  determined 
by  the  amount  of  oxygen  entering  the  porous  medium  at  x  =  0 . 
This  is  converted  to  heat  generation  by  using  the  stoichio- 
metric ratio,  f„,  and  the  heat  of  combustion,  AH_.   Thus  in 
equation  III.  42,  the  planar  heat  source  becomes  puf_AH„(})^. 
Ozisik  [40]  presents  a  general  discussion  and  references  for 
the  moving  heat  source  problem. 

Transition  from  the  internal  combustion  problem  to  the 
surface  recession  problem  occurs  when  the  porosity  at  x  =  0 
nears  a  value  of  1.   The  computer  program,  implementing  the 
analysis,  was  written  so  that  the  change  from  an  internal 
combustion  problem  to  a  surface  recession  problem  occurs 
automatically. 

H.   BOUNDARY  CONDITIONS 

Three  sets  of  boundary  conditions  may  be  used  for  equations 
III. 32,  III. 38,  and  11.39.   Each  set  of  boundary  conditions 
approximates  a  physical  situation  (e.g.,  thermal  flow  reactor 
or  Fontenot's  experiments  [1]).   The  boundary  conditions 
were  the  nearest  approximations  that  could  be  made  and  remain 
within  the  limitations  of  a  one  dimensional  model.   The  first 
and  second  set  of  boundary  conditions  are  typical  of  thermal 
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flow  reactors  with  and  without  radiation  from  the  boundary- 
surf  aces.   These  are 


8T 
c 


=   0 


X  =  0      (III. 45) 


lie 

3x 


=   0 


X  =  L      (III. 46) 


k  -^     =      p  c  u(T   -  T  ) 
a  8x      a  a    a    °° 


X  =  0      (III. 47) 


3T 
a 

3x 


=   0 


X  =  L      (III. 48) 


"e  11=   ^(*  -  *«' 


X  =  0      (III. 49) 


I*   =   0 

3x 


X  =  L      (III. 50) 


and  with  radiation. 


^"^c         "4    ^^4 
(k  +k  )  -^  =   ae  (T   -  T  ) 
e   r   9x         c    <» 


x  =  0      (III. 51) 


ST 


(k  +k  )  ^  =  -aeii^   -   T^) 
e  r   dx  c    «> 


X  =  L      (III. 52) 


ST 

k   -^  =   p  c  u(T   -  T  ) 
a  3x       a  a    a    <» 


X  =  0      (III. 53) 


3T 
a 

3x 


=   0 


X  =  L      (III. 54) 


p   |1  =   u((l)  -  d)  ) 
e  3x       ^^    ^oo' 


X  =  0      (III. 55) 
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1^  =   0  X  =  L      (III. 56) 

dX 


Expressions  III. 4 5  and  III. 46  represent  insulated  boundaries 
on  the  porous  solid  (i.e.,  no  heat  loss  from  the  porous  solid 
to  the  environment).   Expressions  III. 51  and  III. 52  provide 
for  radiation  heat  transfer  from  the  porous  solid  boundaries 
to  the  environment.   The  boundary  conditions  for  the  air 
temperature  and  oxygen  concentration,  expressions  III. 4 7 
through  III. 49  and  I I I. 5 2  through  I I I. 55,  are  the  Danckwerts ' 
boundary  conditions  for  flow  reactors.   A  convincing  discussion 
of  the  Danckwerts'  boundary  conditions  applied  to  mass  diffu- 
sion equations  for  porous  media  is  given  by  Bischoff  [41] . 
A  brief  summary  of  the  discussion  was  presented  by  Vatikiotis 
[2].   An  analogy  to  Bischoff ' s  approach  for  the  air  heat  trans- 
fer equation  is  presented  in  Appendix  B. 

The  third  set  of  boundary  conditions  approximates  the  situa- 
tion for  Fontenot's  experiments  [1].   In  the  experiments,  plate 
specimens  were  mounted  into  the  wall  of  a  wind  tunnel.   The 
interior  surface  of  the  plate  (i.e.,  at  x  =  L)  was  exposed  to 
the  wind  tunnel  air  flow,  and  the  external  plate  surface  (i.e., 
at  X  =  0)  was  exposed  to  the  outside  environment.   The  expres- 
sion for  calculating  the  pressure  differential  across  the 
porous  medium  as  a  result  of  the  external  flow  over  the  surface 
at  X  =  L  is  formulated  in  Appendix  B.   It  must  be  emphasized 
that  the  purpose  for  modelling  Fontenot's  experiments  was  to 
analyze  the  combustion  behavior  within  the  plate  specimens 
using  the  one  dimensional  model.   However,  the  external  flow 
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and  plate  surfaces  involved  convection  heat  transfer  requiring 
two  dimensional  analysis.   Noting  this,  the  following  approach 
was  adopted  to  provide  for  convection  heat  transfer  at  the 
boundaries  of  the  one  dimensional  model  in  a  qualitative 
sense.   The  boundary  conditions  are. 


8T 


(k  +k  )^r-£  =  h^  (T  -T  )+ae{T^-T^) 
e   r  8x     1   c   0°     '  c   « 


X  =  0 


(III. 57) 


9T 


(k  +k  )^r-^  =   -h^(T  -T  )"ae(T^-T^) 
e   r  9x      2   c   «     '  c   «> 


X  =  L 


(III. 58) 


T   =   T 
a      <» 


X  =  0 


(III. 59) 


aT 

a 

3x 


=   0 


X  =  L 


(III. 60) 


(J)   =  .  cj,^ 


X  =  0 


(III. 61) 


|1   =   0 
3x 


X  =  L 


(III. 62) 


As  stated  above,  the  values  of  h,  and  h^  in  the  convection 
terms  of  expressions  III. 57  and  III. 58  are  difficult  to  inter- 
pret for  a  one  dimensional  model.   The  difficulty  arises  from 
considering  the  boundary  layers  on  the  external  surfaces  of 
the  porous  medium  (e.g.,  considering  the  porous  medium  as  a 
flat  porous  plate,  the  boundary  layer  resulting  from  natural 
convection  on  the  x  =  0  surface  increases  in  thickness  in  the 
vertical  direction,  and  for  forced  convection  at  the  x  =  L 
surface,  the  boundary  layer  thickness  increases  in  the  flow 
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direction) .   For  this  investigation,  h,  was  approximated 
by. 


k  u  p 
h,   =   (-^-T7— ^)  (III. 63) 

^    x=0 


This  expression  is  based  on  the  analytical  results  presented 
by  Merkin  [42].   Expression  III. 63  is  obtained  by  considering 
natural  convection  heat  transfer  from  a  vertical  flat  plate 
with  "suction".   From  Merkin ' s  results,  the  boundary  layer 
thickness  becomes  constant  at  some  distance  in  the  vertical 
direction,  and  thus,  the  convection  heat  transfer  may  also  be 
considered  constant.   The  magnitude  of  the  convection  heat 
transfer  coefficient,  h^,  varies  in  a  direction  parallel  to 
the  external  flow,  U   at  the  x  =  L  surface.   In  addition,  h- 
depends  on  the  efflux  of  the  gas  at  the  surface.   To  simplify 
the  analysis,  h^   was  approximated  by  the  relations  for  a  smooth 
flat  plate  given  by  Holman  [35]  as. 


h^      =   .664  -p  Pr-^^-^  Re-^/^  (III. 64) 


for  laminar  flow,  where  L  is  the  distance  from  the  boundary 
layer  initiation,  and. 


h^      =      -p  Pr-^/^(.037  Re*^  -  850)  (III. 65) 


for  turbulent  flow,  where  Pr  is  the  Prandtl  number.   The 
Reynolds  number  here  is  based  on  the  external  flow  velocity 
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over  the  surface  at  x  =  L.   Kays  [43]  provides  an  alternate 
scheme  for  treating  boundary  layers  with  "blowing  and  suction", 
The  air  temperature  and  oxygen  concentration  boundary  condi- 
tions, III. 59  through  III. 62  were  also  selected  because  of 
the  difficulties  in  treating  external  surface  boundary  layers 
in  a  one  dimensional  model .   Alternatives  to  the  essential 
boundary  conditions  of  the  air  and  oxygen  are  the  Danckwerts ' 
conditions.   However,  using  the  Danckwerts'  conditions  with 
convection  boundary  conditions  poses  the  problem  of  approxi- 
mating a  reasonable  value  for  u  at  x  =  0  and  x  =  L  (i.e., 
boundary  layers  affect  the  component  of  velocity,  u,  normal 
to  an  external  surface).   To  overcome  these  difficulties,  it 
seemed  prudent  not  to  have  pore  velocity  and  the  convection 
heat  transfer  coefficients,  h,  and  h^ ,  appear  in  the  same  set 
of  boundary  conditions.   However,  the  difficulty  still  re- 
mains of  specifying  a  distance  from  the  origin  of  the  boundary 
layer  for  expressions  III. 64  and  III. 65.   In  the  initial  work 
by  Vatikiotis  [2] ,  an  arbitrary  distance  of  1  foot  was  used 
to  compare  results. 

I.   INITIAL  CONDITIONS 

To  initiate  combustion,  the  porous  medium  must  be  brought 
to  a  temperature  at  which  the  heat  generation  rate  is  suffi- 
ciently high.   Experimentally,  there  are  many  techniques  for 
doing  this.   However,  it  is  difficult  to  measure  the  actual 
temperature  and  oxygen  concentration  during  the  experiments. 
This  difficulty  is  carried  over  to  identifying  a  reasonable 
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set  of  initial  conditions.   Noting  this,  the  model  was  de- 
veloped so  a  problem  can  be  started  at  ambient  conditions 
with  a  heat  flux  applied  to  the  carbon  at  a  boundary.   The 
heat  flux  is  mathematically  treated  as. 


'l-P"^  +  "^r*  ^  =  -%t  (III. 66) 


This  approach  of  starting  the  problem  simplifies  the  task  of 
determining  a  reasonable  set  of  initial  conditions.   In  addi- 
tion, any  arbitrary  set  of  initial  conditions  can  be  specified, 
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IV.   IMPLEMENTATION  OF  NUMERICAL  METHODS 

A.   SOLUTION  ALGORITHM 

The  field  equations  and  auxiliary  equations  which  define 

the  problem  were  presented  in  the  previous  sections.   The 

following  is  the  scheme  used  to  solve  for  the  dependent 

variables  T  ,  T  ,    <^ ,    u,  P,  L,  and  p  : 
c   a   ^ '   '   '         a 

1.  Starting  with  the  Dupuit-Forcheimer  expression  III. 10, 
the  filter  velocity,  Q,  is  represented  as  a  function 
of  the  pore  velocity,  u. 

2.  Substitute  for  Q  from  the  previous  step  into  Darcy's 
law,  equation  A.l. 

3.  The  expression  formed  in  step  2  is  solved  for  the  pore 
velocity,  u,  as  a  function  of  dP/dx.   This  procedure 
yields  equation  A. 2. 

4.  Equation  A. 2  is  substituted  into  the  continuity  equation 

A.  3  resulting  in  the  2nd  order  partial  differential 

equation  A. 4  with  the  pressure,  P,  and  the  air  density, 

p  ,  as  the  dependent  variables. 
a 

5.  Expanding  equation  A. 4  yields  equation  A. 5.   Equation 
A.  5  is  one  of  4  equations  which  were  integrated  to 
produce  the  problem  solution.   By  taking  the  air  den- 
sity, p  ,  as  a  constant  over  a  time  step  of  integration, 

a 

equation  A. 5  becomes  an  ordinary  differential  equation 
with  respect  to  x.   For  equation  A. 5,  the  shooting  method 
(with  Euler's  formula  for  the  integration  algorithm) 
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was  used  to  solve  for  the  pressure,  P (x) ,  and  the 

pressure  gradient,  dP(x)/dx  at  each  time  step.   The 

air  density,  p  ,  was  updated  at  each  time  step  using 
a 

the  ideal  gas  law,  equation  B.41.   The  rate  of  change 
of  the  air  density  was  neglected  during  the  first  two 
time  steps  of  integration.   Reasons  for  neglecting  the 
time  rate  of  change  of  the  air  density  are  given  in 
Appendix  C.2. 

6.  Having  solved  equation  A. 5  for  the  pressure  gradient, 
dP/dx,  the  pore  velocity  can  be  obtained  from  equation 
A. 2. 

7.  In  solving  the  air  heat  transfer  equation  III. 38  and 

the  oxygen  molecule  transfer  equation  III. 39,  the  pore 

velocity,  u,  and  the  air  density,  p  ,  are  obtained  from 

a 

equations  A. 2  and  B.41,  respectively.   Equation  III. 38 
and  III. 39,  along  with  the  carbon  particle  equation 
III.  32,  were  transformed  to  ordinary  differential  equa- 
tions in  time  by  a  Galerkin  FEM  formulation. 

8.  The  ordinary  differential  equations  from  step  7  are 
given  by  equations  C.17  through  C.19.   The  time  inte- 
gration was  performed  by  a  version  of  Gear's  method 
developed  by  Franke  [44] . 

9.  All  thermophysical  properties  (i  .e .  ,  k  ,  k  ,  k  ,  c  , 

c  ,  y)  being  functions  of  temperature  are  continuously 
updated  during  the  transient  analysis. 
10.   As  carbon  is  consumed  during  combustion,  the  porous 

medium  properties  (i.e.,  p,  m,  5,  d)  change  with  time 
and  are  updated  continuously  as  well. 
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12.   No  attempt  was  made  to  match  the  orders  of  convergence 
of  the  integration  algorithms. 

Implementing  the  above  solution  algorithm  resulted  in  a 
computer  program  for  the  solution  of  the  problem.   Several 
options  were  incorporated  in  the  computer  program  allowing 
the  user  flexibility  in  the  analysis.   Boundary  conditions 
simulating  fixed  air  inlet  temperature,  and  convection  and 
radiation  heat  transfer  from  the  exterior  surfaces  can  be  se- 
lected.  The  problem  can  be  initiated  in  two  ways;  an  arbitrary 
set  of  initial  conditions  can  be  specified,  or  a  heat  flux, 
shown  by  expression  III. 66,  can  be  specified  with  the  porous 
medium  at  ambient  conditions.   In  addition,  the  porosity  at 
X  =  0  causing  transition  to  the  surface  recession  phase  can 
be  specified.   Other  options  permit  investigating  materials 
other  than  carbon  (e.g.,  boron),  studying  the  effects  of  non- 
constant  distributions  of  porous  media  properties. 

The  computer  program  includes  a  main  program  and  24  sub- 
routines (7  of  the  subroutines  comprise  Franke ' s  [44]  inte- 
gration routine) .   The  function  and  description  of  each  of 
the  subroutines  are  found  in  the  computer  program  listing. 
Figure  IV. 1  shows  the  flow  of  the  program,  and  briefly  des- 
cribes the  process  at  each  step.   Each  block  represents  a 
subroutine.   The  configuration  of  the  blocks  (i.e.,  blocks 
within  blocks)  in  Figure  IV. 1  show  the  six  levels  of  the 
computer  program.   In  addition,  other  processes  associated 
with  the  integration  routine  occur  inside  the  block  shown  by 
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Main  program. 


Read 

in  p 

arameters, 

ini  t 

i  al  i 

se 

■fl 

ags 

and 

coun 

ters; 

Calculate 

ini 

tial 

prop 

ert 

ies 

and 

pore 

vel oci ty. 

Construct  FEM  nodal  point/element 
correspondence  array. 


Print  current  response  variables 
and  properties. 


Check  for  transition 
to  sur-face  recession. 

TRUE 

Reset  -flags,  counters  delete 

<t>   (cal  led  once)  . 

1 
FALSE 

1 


Call  integration  routine  -For  T  ,  T  ,  d)  ,  (  C±)  deleted 
during  sur-face  recession). 

^Construct  O.D.E.'s. 

Cal 

culate  properties. 

Test  -for  oxygen  instabiltiy. 

I-f  sur-face  recession,  update  thickness,  L. 

Calculate  pore  velocity  and  pressure  gradient. 

I-f  sur-face  recession,  trans-form  reaction  rate 
expression  to  sur-face  -flux. 

Set  -flag  -for  sur-face  recession  i-f  p(n=l)  =  pma; 

... 

Zero-out  system  matrices. 

Generate  FEM  operators. 

Form  mass  and  system  matrices. 

Adjust  matrices  -for  boundary  conditions. 

Adjust  matrices  -for  oxygen  instability. 

jPer 

■■form  integration. 

FALSE — Test    -flags    for    program    termination. 


TRUE 


Print  last  values  of  response  variables 
and  properties. 


STOP 


FIGURE  IV. 1  -  Flow 
of  computer  program, 


55 


the  dashed  line.   These  processes  are  discussed  in  detail  in 
Franke's  report. 

B.   STIFFNESS  CONSIDERATIONS 

Initially,  the  ordinary  differential  equations  resulting 
from  the  finite  element  formulation  were  integrated  explicitly 

using  a  sixth  order  Runge-Kutta  method  (IMSL  subroutine  DVERK) . 

-4 
The  largest  obtainable  time  step  was  approximately  10 

seconds.   A  larger  time  step  caused  the  integration  to  become 
unstable.   Gear  [45]  described  this  behavior  as  typical  of 
"stiff"  differential  equations,  and  suggested  that  implicit 
integration  methods  be  used  to  improve  the  efficiency  of  the 
integration.   Gear's  method,  presented  by  Brown  and  Gear  [4  6] 
and  modified  by  Franke  [44],  was  adopted  as  the  integration 
method.   As  a  result,  time  steps  of  several  seconds  were  ob- 
tained for  some  problems  (i.e.,  for  surface  recession  and  when 
extinction  occurred) .   In  addition,  moving  the  exponential 
reaction  rate  from  the  excitation  vector  to  the  stiffness 
matrix  (discussed  in  Appendix  C.l)  improved  the  integration 
efficiency.   The  improved  efficiency  was  attributed  to  a  better 
approximation  of  the  Jacobian  matrix  in  Gear's  method. 

Bui  [47]  and  Burka  [48]  discuss  the  difficulties  of  solving 
differential  equations  exhibiting  "stiffness"  and  propose 
alternatives  to  Gear's  method.   For  a  system  exhibiting  "stiff- 
ness", small  time  steps  are  needed  for  the  rapidly  responding 
terms  while  the  integration  must  be  continued  for  a  long  period 
to  account  for  the  slowly  responding  terms.   A  measure  of 


56 


"stiffness"  for  a  system  of  differential  equations  is  the 
ratio  of  the  real  part  of  the  maximiam  eigenvalue  of  the  Jacob- 
ian  matrix  to  the  real  part  of  the  minimum  eigenvalue.   This 
is  called  the  "local  stiffness  ratio",  and  is  represented 
by. 

Max      I  A .  I 
i 

5  =  ^"?;;;"^" — n—r  (iv.i) 

Mm      1^.1 
i=l , ...  ,n 

where  n  is  the  number  of  equations.   A  system  of  equations  can 
be  considered  "stiff"  if  S  has  a  value  greater  than  10,  and 
the  real  parts  of  the  eigenvalues  are  negative. 

Lastly,  most  numerical  integration  methods  are  bound  by 
a  maximum  time  step  for  which  the  solution  remains  stable. 
This  is  defined  by. 


At-Aj_|   <   K  (IV. 2) 


where  At  is  the  time  step.  A-  is  the  real  part  of  the  maximum 
eigenvalue  of  the  Jacogian  matrix,  and  K  is  a  constant  depend- 
ing on  the  integration  method  (usually  1  <  K  <  10) .   For 
example,  Bui  [47]  states  that  Euler's  method  requires 
|AtA.  j  <  1.   Therefore,  for  a  maximum  Jacobian  matrix  eigen- 
value (representing  fastest  response)  of  50  (seconds)"  ,  the 
largest  timestep  for  the  integration  is  .02  seconds.   Exceed- 
ing this  value  results  in  the  integration  becoming  unstable. 
This  explains  the  maximum  timestep  of  lO"   seconds  for  the 
sixth  order  Runge-Kutta  method  as  previously  discussed. 
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C.   TREATMENT  OF  NUMERICAL  DIFFICULTIES  ASSOCIATED  WITH 
OXYGEN  CONCENTRATION 

In  problems  of  combustion  or  kinetics,  where  dependent 
variables  go  to  zero,  numerical  difficulties  arise.   In  his 
discussion,  Frank-Kamenetskii  [12]  points  out  that  for  frac- 
tional order  reactions  (expression  III. 16),  both  the  oxygen 
concentration  and  concentration  gradient  can  become  equal  to 
zero  within  the  porous  medium.   This  requirement  is  difficult 
to  satisfy  with  approximate  solution  methods.   In  the  initial 
effort  (Vatikiotis  [2]),  run  time  errors  resulted  when  the 
oxygen  concentration  became  slightly  negative.   This  was 
caused  by  attempting  to  take  the  logarithm  of  a  negative  con- 
centration in  the  fractional  order  term  in  expression  III. 16. 
One  approach  for  correcting  this  problem  was  to  use  the  abso- 
lute value  of  the  oxygen  concentration  in  the  reaction  rate 
term.   This  resulted  in  the  concentration  becoming  increasingly 
negative.   Another  approach  was  to  check  the  values  of  oxygen 
concentration  at  each  integration  interval  and  setting  the 
concentration  to  zero  for  the  negative  values.   At  the  next 
integration,  the  reaction  rate  was  zero  at  the  locations  of 
zero  concentration.   Without  consumption  the  oxygen  concentra- 
tion at  the  nodal  point  previously  set  to  zero  becomes  positive 
and  the  oscillation  repeats  itself.   It  was  observed  that  de- 
creasing the  integration  time  step  and  the  length  of  the  ele- 
ments (i.e.,  increasing  the  number  of  nodal  points)  minimized 
the  oscillations.   Since  CPU  time  increased  for  a  given  problem, 
the  approach  was  not  totally  satisfactory.   However,  this 
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method  was  adopted  for  obtaining  the  results  presented  by 
Vatikiotis  [2] . 

An  attempt  to  resolve  the  above  numerical  difficulty  for 
the  present  model  was  by  implementing  the  Moving  Finite  Element 
method  of  Gelinas,  Doss,  and  Miller  [49] .   In  the  method, 
the  nodal  points  defining  the  elements  converge  to  regions 
where  they  are  needed  to  minimize  the  error  in  the  numerical 
solution.   As  a  result,  the  size  of  the  elements  in  a  region 
of  large  reaction  rate  become  small,  and  move  with  this  region 
in  time.   The  motivation  for  implementing  the  Moving  Finite 
Element  method  was  Frank-Kamenetskii ' s  [12]  discussion  of 
fractional  order  reactions  as  stated  above.   Also,  it  was  ob- 
served that  the  oxygen  concentration  approached  zero  in  a 
small  interval  within  the  porous  medium.   This  behavior  resulted 
in  steep  concentration  gradients.   It  was  thought  that  a  num- 
ber of  small  elements  in  this  region  would  improve  the  accuracy 
of  the  solution.   The  results  obtained  from  the  Moving  Finite 
Elements  showed  the  method  was  not  effective  in  resolving  the 
numerical  difficulties  discussed  above.   The  method  was  subse- 
quently abandoned  because  of  the  additional  expense  associated 
with  solving  the  differential  equation  governing  the  nodal 
point  locations. 

The  method  adopted  in  the  present  model  for  minimizing 
the  difficulties  with  the  oxygen  concentration  is  as  follows. 
The  time  derivative  of  the  oxygen  concentration  at  a  nodal 
point  is  set  to  zero  when  the  concentration  reaches  a  speci- 
fied  positive  minimum  value  (less  than  10    Ibm/ft  ) .   At 
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the  next  integration  interval,  the  time  derivative  remains 
zero  if  the  concentration  at  the  adjacent  nodal  point  in  the 
negative  direction  (i.e.,  upstream)  has  decreased.   If  the 
concentration  at  the  adjacent  nodal  point  has  increased,  then 
the  time  derivative  is  released.   This  method,  in  conjunction 
with  Franke ' s  [44]  integration  routine,  appears  to  act  in 
an  iterative  manner  continually  improving  the  solution.   The 
values  of  concentration  in  the  regions  where  the  oxygen  has 

been  totally  consumed  are  typically  on  the  order  of  10 

3 
Ibm/ft  .   Although  they  do  not  give  the  details,  Scheisser 

and  Stein  [50]  have  used  the  above  method  in  their  worlc  of 

simulating  coal  conversion. 

D.   AN  ALTERNATE  SOLUTION  STRATEGY 

The  experience  and  difficulties  of  implementing  the  numeri- 
cal methods,  as  discussed  above,  have  provided  insight  for  an 
alternate  strategy  for  solving  the  combustion  and  heat  trans- 
fer problem.   The  alternate  strategy  consists  of  formulating 
the  problem  with  two  moving  boundaries.   One  moving  boundary, 
as  in  the  present  formulation,  would  account  for  surface 
recession  (i.e.,  the  porous  medium  thickness  changing  with 
time) .   The  second  moving  boundary  would  account  for  the 
oxygen  concentration  going  to  zero  within  the  porous  medium. 
This  second  boundary  would  be  located  where  the  concentration 
becomes  zero,  between  x  =  0  and  the  first  moving  boundary. 
The  advantages  of  the  alternate  solution  strategy  are  (1) 
eliminating  the  numerical  difficulties  associated  with  the 
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oxygen  concentration  (discussed  in  the  previous  section) , 
(2)  relaxing  the  planar  heat  source  assumption  for  the  surface 
recession  phase  (discussed  in  Section  III.G),  and  (3)  possibly 
improving  the  numerical  "stiffness"  characteristics  of  the 
integration  (discussed  in  Section  IV. B) . 
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V.   RESULTS  AND  OBSERVATIONS 

A.  BASIS  FOR  THE  RESULTS 

The  results  presented  in  the  following  sections  are  based 
on  the  idealized  porous  medium  described  in  Section  III. A  and 
Figure  III.l.   For  the  idealized  porous  medium,  porosity  and 
permeability  were  explicit  functions  of  pore  diameter  and  unit 
cell  thickness.   Other  porous  media  will  have  different  geo- 
metric relationships  for  describing  the  properties.   In  addi- 
tion, for  naturally  occurring  porous  media,  experimental  methods 
would  be  employed  to  determine  the  geometric  properties. 
Noting  this,  the  results  of  this  investigation  are  intended 
to  describe,  in  a  qualitative  sense,  the  behavior  of  a  carbon 
porous  medium  during  combustion.   As  will  be  shown,  the  tempera- 
tures observed  during  combustion  are  highly  dependent  on  the 
geometric  properties  and  the  pore  velocity.   Although,  the 
temperatiires  are  typical  of  those  reported  for  carbon  combustion 
experiments,  it  would  not  be  reasonable  to  perform  a  quanti- 
tative comparison  between  the  results  obtained  by  the  model 
and  those  obtained  experimentally  unless  the  porous  medium 
properties  were  the  same  for  both.   This  point  will  be  demon- 
strated in  the  next  section. 

B.  COMPARISON  OF  COMBUSTION  MODEL  RESULTS  TO  EXPERIMENTAL 
RESULTS 

The  behavior  of  the  temperature  and  oxygen  concentration 

of  the  combustion  and  heat  transfer  model  is  in  fair  agreement 
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with  the  behavior  reported  in  the  experimental  investigations 
of  Vulis  [13] ,  Koizumi  [37] ,  Fontenot  [51] ,  Spalding  [52] ,  and 
Kolodstev  [53] .   A  quantitative  comparison  of  the  results  can 
not  be  made  for  two  reasons.   First,  there  have  been  few 
experimental  investigations  which  present  data  for  transient 
behavior  of  combustion  in  porous  media.   Most  of  the  analyses 
present  quasi-steady  data  (quasi  in  the  sense  that  for  some 
period  of  time,  there  is  little  change  in  the  response  varia- 
bles) .   When  transient  behavior  is  reported,  it  is  usually 
presented  in  a  narrative  form.   Secondly,  the  reported  des- 
criptions of  the  porous  media  used  in  the  experiments  are  not 
sufficient.   The  properties  of  the  carbon  used  (e.g.,  thermal 
conductivity,  density,  Arrhenius  rate  law,  etc.)  and  the 
porosity  and  permeability  relationships  are  needed  to  accurately 
simulate  the  experiment.   The  profiles  of  the  quasi-steady 
temperatures  obtained  by  the  combustion  model  and  those  ob- 
tained by  Kolodstev* s  experiments  will  be  shown  together.   The 
purpose  of  this  is  to  demonstrate  that  although  Kolodstev' s 
porous  medium  could  only  be  roughly  approximated,  the  tempera- 
ture profile  obtained  by  the  combustion  model  was  similar  in 
shape  and  magnitude  to  that  obtained  experimentally. 

In  the  experimental  investigations  of  Vulis  [13] ,  Fontenot 
[51],  and  Spalding  [52] ,  it  was  observed  that  for  a  given 
temperature,  sustained  combustion  occurred  for  a  specific  range 
of  velocities.   Velocities  outside  the  range,  either  greater 
or  lower,  would  cause  the  combustion  to  extinguish  and  the 


63 


porous  medium  to  cool  to  ambient  temperature.   This  phenomenon 
was  also  observed  for  the  combustion  model.   The  explanation 
of  this  behavior,  discussed  in  Section  V.G,  is  summarized  as 
follows.   As  pore  velocity  increases,  convection  heat  transfer 
and  oxygen  supply  for  combustion  both  increase.   Similarly, 
a  decrease  in  pore  velocity  will  cause  a  decrease  in  convec- 
tion heat  transfer  and  oxygen  supply.   For  a  given  set  of 
conditions,  pore  velocities  above  the  range  will  result  in 
the  heat  transfer  or  heat  loss  overcoming  the  heat  generation. 
Also,  for  pore  velocities  below  the  range,  the  oxygen  supply 
for  combustion  is  decreased  to  the  extent  that  the  heat  genera- 
tion is  overcome  by  the  heat  loss.   In  other  words,  the  higher 
pore  velocities  tend  to  "blowout"  the  reaction  (e.g.,  blowing 
out  a  candle) ,  and  the  lower  pore  velocities  tend  to  "starve" 
the  reaction  (e.g.,  placing  a  container  over  a  candle). 
Important  at  lower  pore  velocities  are  the  effects  of  heat 
transfer  from  the  external  surface  of  the  porous  solid  (i.e., 
at  the  boundaries) .   The  experiments  performed  by  Vulis  [13] , 
Fontenot  [51] ,  and  Spalding  [52]  allowed  heat  transfer  from 
the  surfaces  by  either  convection  or  radiation.   The  effects 
of  heat  transfer  at  the  boundaries  on  combustion  are  discussed 
in  Section  V.G. 

Kolodstev  [53]  investigated  combustion  gas  dynamics  in  a 
porous  medium  comprised  of  spherically  shaped  carbon  particles. 
This  is  the  same  type  of  carbon  used  by  Parker  and  Hottel  [16] . 
In  his  results,  Kolodstev  presents  the  quasi-steady  air  tempera- 
ture as  a  function  of  the  depth  of  the  porous  medium  (particle 
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temperature  was  not  reported).   Figure  V.l  shows  Kolodstev's 
results  plotted  with  the  results  obtained  by  the  combustion 
model.   Kolodstev's  experiment  could  not  be  simulated  accu- 
rately.  As  stated  previously,  the  description  of  Kolodstev's 
porous  medium  was  not  sufficient.   The  description  of  the 
porous  medium  and  the  ambient  conditions  used  to  simulate  the 
experiment  is  shown  in  Table  V.l.   Also  shown  are  the  known 
parameters  given  by  Kolodstev.   In  Figure  V.l,  the  temperature 
for  the  first  inch  of  the  7 . 5  inches  of  the  porous  medium  are 
shown  (temper a tiires  at  greater  depths  were  not  reported)  .   As 
can  be  seen,  there  is  fair  agreement  in  the  shapes  of  the 
profiles.   In  addition,  the  magnitude  of  the  temperatures  pro- 
duced by  the  combustion  model  are  similar  to  those  obtained 
in  the  experiment. 

The  combustion  model  results  showed  that  once  the  reaction 
moved  to  the  air  inlet  surface  of  the  porous  medium,  the  oxy- 
gen concentration  penetrated  the  porous  medium  by  1  or  2 
particle  diameters.   This  observation  was  also  reported  in 
the  experimental  results  of  Koizumi  [37]  and  Kolodstev  [53] . 
The  response  of  the  oxygen  concentration,  in  general,  will  be 
discussed  in  the  next  section. 

C.   EXAMPLE  RESULTS 

The  intent  here  is  simply  to  demonstrate  the  operation  of 
the  computer  program.   The  behavior  of  a  single  porous  medium 
subjected  to  several  initial  conditions  is  presented.   The 
description  of  the  porous  medium  and  the  ambient  conditions 
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TABLE  V.  1 Geometry  o-f  porous  medium  and  ambient  conditions 

for  simulating  Kolodstev's  C53]  experiment. 


•  Particle  Shape:  spherical 

•  Particle  Diameter  (  d  )  :  .126  in 
Unit  Cell  Thickness  (Q):  .126  in 

•  Spatial  Thickness  o-f  Porous  Medium  (L  ^i  7.5  in 

Porosity  ( p  )  .476 

-A     2 

Permeability  (  m  )  .103x10^  -Ft 

Bulk  Thermal  Conductivity  o-f  Carbon  (  k^  )  :  86.0  Btu/-ft-hr — F 

Bulk  Speci-fic  Heat  of  Carbon  (  C^ )  :  .231  Btu/lbm-F 

Bulk  Density  o-f  Carbon  Cp  ):  70.3  lbm/-ft 

Thermal  Emissivity  o-f  Particles  (€):  .9 

Ambient  Temperature  (  X^) :  80.0  deg-F 

•  Ambient  Pressure  (F^):  14.7  psi 

•  Pressure  Di-f -f  erenti  al  (aP)  :  -.29  psi 

•  Ambient  Oxygen  Concentration  (ct>.)  :  .0172  Ibm/ft 


•  Known  parameter  -for  Kolodstev's  experiment. 
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used  for  the  following  problems  are  shown  in  Table  V.2.   The 
boundary  conditions  for  the  examples  in  this  section  were 
insulated  boundaries  on  the  porous  solid  temperature  (ex- 
pressions III. 45  and  III. 46),  and  the  Dankwerts  conditions  on 
the  air  temperature  and  oxygen  concentration  (expressions 
III. 47  through  III. 50). 

As  the  first  example  (referred  to  as  example  A) ,    a  heat 

4       2 
flux  of  3.0  X  10   Btu/ft  -hr  was  applied  at  x  =  0  for  15  seconds 

The  porous  medium  was  initially  at  a  constant  temperature  of 

80  degrees  Fahrenheit  and  a  constant  oxygen  concentration  of 

3 
.0172  Ibm/ft   (ambient  conditions) .   The  results  of  this  prob- 
lem are  shown  in  Figures  V.2-V.4.   Figure  V.2  shows  the  tem- 
perature increasing  after  the  heat  flux  was  removed.   At  25 
seconds,  the  porosity  at  x  =  0  reached  .95  and  the  problem 
transitioned  to  the  surface  recession  phase.   The  values  of 
thickness  for  the  times  shown  in  Figure  V.2  are  .25  inches  to 
25  seconds,  .226  inches  at  72  seconds,  .133  inches  at  357 
seconds,  and  .0238  inches  at  517  seconds.   The  computer  pro- 
gram was  written  such  that  the  problem  terminates  when  the 
thickness  becomes  10  percent  or  less  of  the  initial  value. 
The  air  temperature  shown  by  the  dashed  line  in  Figure  V.2 
follows  the  carbon  temperature  with  the  exception  of  the  air 
inlet  region  at  x  =  0 .   Figure  V.3  shows  the  oxygen  behavior 

for  this  problem.   At  25  seconds,  the  oxygen  concentration  at 

-4       3 
X  =  0  was  3.6  X 10    Ibm/ft   and  the  penetration  of  the  oxygen 

was  to  x/L  =  .01.   This  supports  the  assumption  of  using  a 

planar  source  for  the  heat  generation  and  deleting  the  oxygen 
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TABLE  V.2 Geometry  o-f  porous  medium  and  ambient  conditioni 

■for  the  example  problems  in  Section  V.C. 


Particle  Shape:  spherical 

Particle  Diameter  (d):  .005  in 

Unit  Cell  Thickness  (  D):  .005  in 

Spatial  Thickness  o-f  Porous  Medium  (  L):  .25  in 

Porosity  ( p  )  . 476 

Permeability  (  m  )  .162x10^  -ft 

Bulk  Thermal  Conductivity  o-f  Carbon  (  k^.)  :     86.0  Btu/-ft-hr — F 

Bulk  Speci-fic  Heat  o-f  Carbon  (  C^  )  :              .231  Btu/lbm-F 

Bulk  Density  of  Carbon  (p^  ^  *  70.3  lbm/-ft 

Thermal  Emissivity  o-f  Particles  (€):  .9 

Ambient  Temperature  (T  >!  SO.O  deg-F 

Ambient  Pressure  (  f^)  :  14.7  psi 

Pressure  Di-f -f  erenti  al  (Ap)  :  -.35  psi 


Ambient  Oxygen  Concentration  <  <t>^)  :  .0172  lbm/-ft 
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concentration  equation  during  the  surface  recession  phase. 
Figure  V.4  shows  the  decrease  in  the  thickness  of  the  porous 
medium  with  time.   The  nondimensional  parameter,  X,    given  by 
expression  III. 31,  had  an  average  value  of  .63  at  the  start  of 
the  problem.   This  value  of  X   is  above  the  minimum  value  of 
.342  proposed  by  Green  and  Perry  [20]  when  assuming  equal 
carbon  and  air  temperatures.   This  suggests  that  a  single  heat 
transfer  equation  would  have  been  sufficient  to  describe  the 
temperature  response  for  this  problem.   However,  this  is  not 
valid  for  all  cases.   For  example,  when  d,  D,  and  L  are 
doubled,  that  is,  .01  inches,  .01  inches,  and  .5  inches, 
respectively,  A  equals  .22,  and  for  d,  D,  and  L  equal  to  .025 
inches,  .0  25  inches  and  1.25  inches,  respectively,  X   equals 
.056.   Hence,  for  the  general  case,  the  carbon  and  the  air 
temperatures  should  be  treated  as  separate  response  variables. 

As  a  second  example  (referred  to  as  example  B) ,  a  heat 

4       2 
flux  of  3.0  X  10   Btu/ft  -hr  was  applied  at  x  =  0  for  14  seconds 

The  porous  medium  was  initially  at  the  same  ambient  tempera- 
ture and  oxygen  concentration  as  in  the  previous  example. 
Figures  V.5  and  V.6  show  the  results  of  this  problem.   The 
porous  medium  cooled  to  ambient  conditions  in  approximately 
41  seconds  after  the  heat  flux  was  removed.   Applying  the 
heat  flux  for  14  seconds  was  not  sufficient  to  produce  the 
conditions  needed  to  sustain  combustion.   However,  increasing 
the  duration  of  the  heat  flux  for  an  additional  second  results 
in  sustained  combustion.   Figures  V.5  and  V.6  show  the  results 
starting  at  12  seconds.   The  behavior  of  the  system  prior  to 
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12  seconds  is  the  same  as  in  Figures  V.2  and  V.3  of  the  pre- 
vious problem. 

As  the  last  example  (referred  to  as  example  C) ,  the  problem 
was  started  with  the  porous  medium  at  constant  temperature 

and  oxygen  concentration  of  1100  degrees  Fahrenheit  and  .006 

3 
Ibm/ft  ,  respectively.   Figure  V.7  shows  the  large  cooling 

effect  of  the  air  entering  the  porous  medium  at  x  =  0  during 

the  early  part  of  the  problem.   As  a  result  of  this  cooling, 

the  heat  generation  began  to  raise  the  temperature  interior 

to  the  porous  medium  before  moving  as  a  wave  to  the  air  inlet 

surface.   The  problem  was  terminated  at  49  seconds  with  a 

porosity  at  x  =  0  of  .96.   This  was  at  the  point  of  transition 

to  the  surface  recession  phase.   Figure  V.8  shows  the  behavior 

of  the  oxygen  for  this  problem.   At  49  seconds,  the  oxygen 

-4       3 
concentration  at  x  =  0  was  4.1 xio    Ibm/ft  and  the  penetra- 
tion of  the  oxygen  was  to  x/L  =  .01.   This  again  supports  the 
assumption  of  treating  the  heat  generation  as  a  planar  source 
and  deleting  the  oxygen  concentration  equation  during  the 
surface  recession  phase. 

Other  data  generated  by  the  model  are  spatial  distributions 
of  the  following  internal  properties:   particle  diameter  (d) , 
porosity  (p) ,  specific  permeability  (m) ,  pressure  (P) ,  pres- 
sure gradient,  pore  velocity  (u) ,  Reynolds  number  (Re) ,  heat 
transfer  coefficient  (h) ,  and  surface  area  per  unit  volume  (z) . 
As  previously  stated  in  the  model  development,  these  system 
properties  are  functions  of  temperature,  and  hence,  change 
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with  time.   In  addition,  since  surface  recession  may  occur, 
the  spatial  thickness  (L)  of  the  porous  medium  is  provided 
as  output  data. 

D.   EFFECT  OF  GEOMETRIC  PARAMETERS 

In  the  initial  effort  by  Vatikiotis  [2] ,    it  was  observed 
that,  for  a  given  set  of  boundary  and  initial  conditions,  the 
internal  geometry  of  a  porous  medium  significantly  influenced 
combustion.   Results  were  reported  describing  the  dependence 
of  combustion  on  particle  size,  and  the  thickness  of  the 
porous  medium.   As  a  refinement  of  that  initial  work,  results 
are  presented  below  which  show  the  influence  of  the  geometric 
parameters  (i.e.,  permeability,  porosity,  porous  medium  thick- 
ness) on  the  minimum  initial  temperature  needed  to  sustain 
combustion.   "Minimum  initial  temperature"  is  defined  here  as 
the  minimiom  constant-value  temperature  distribution  (i.e., 
slope  equal  to  zero)  that  will  sustain  combustion  when  used 
as  an  initial  condition.   A  lower  temperature  will  cause  the 
porous  medium  to  cool  to  ambient  temperature.   As  will  be 
shown  in  Section  V.E,  both  the  shape  and  magnitude  of  the 
temperature  initial  conditions  influences  the  combustion  proc- 
ess.  In  order  to  form  a  basis  for  comparison,  the  analysis 
was  performed  with  uniform  distributions  for  the  initial  con- 
ditions (the  carbon  and  air  having  the  same  initial  tempera- 
ture) .   The  initial  condition  for  the  oxygen  concentration 
(also  uniform)  was  taken  as  the  concentration  found  in  air  at 
the  initial  condition  temperature  and  ambient  pressure.   However, 
any  oxygen  concentration  distribution,  as  an  initial  condition. 
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would  not  have  changed  the  results .   The  boundary  conditions 
used  in  the  analyses  were  the  insulated  boundaries,  expressions 
III. 45  and  III. 46,  for  the  carbon  temperature,  and  the  Danck- 
werts  conditions,  expressions  III. 47  through  III. 50,  for  the 
air  temperature  and  oxygen  concentration.   The  results  of  the 
analysis  were  obtained  by  varying  the  initial  conditions  until 
the  "minimum  initial  temperature"  for  sustaining  combustion 
was  obtained.   A  geometric  parameter  was  then  changed  to  a  new 
value  while  all  other  parameters  remained  the  same.   A  "mini- 
mum initial  temperature"  for  this  problem  was  obtained.   This 
procedure  was  repeated  for  several  values  of  the  same  geometric 
parameter  while  keeping  all  other  parameters  fixed.   The  re- 
sults obtained  showed  the  "minimum  initial  temperature"  which 
sustained  combustion  as  a  function  of  the  respective  parameter. 
The  results  for  each  geometric  parameter  are  discussed 
individually. 

1.   Effects  of  Permeability 

Figure  V.9   shows  the  dependence  of  the  "minimum  initial 
temperature"  (defined  previously  in  this  section)  on  permea- 
bility.  The  description  of  the  porous  medium  and  the  ambient 
conditions  used  in  the  following  analysis  is  given  in  Table 
V.3.   The  range  of  permeability  was  determined  by  the  magni- 
tude of  the  Reynolds  number.   The  largest  permeability 

—  8    2 
(m  =  .4x10    f t  ) ,  at  a  pressure  differential  of  1.5  p.s.i., 

resulted  in  an  average  Reynolds  number  through  the  porous 

medium  of  approximately  5.   It  was  felt  that  extending  the 
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TABLE  V.3 Geometry  o-f  porous  medium  and  ambient  conditions 

■For  the  permeability  analysis  in  Section  V.D. 


Particle  Shape:  spherical 

Particle  Diameter  ( d  ) :  (various) 

Unit  Cell  Thickness  (0^=  (various) 

Spatial  Thickness  o-f  Porous  Medium  (  L  >!  1.0  in 

Porosity  ( p )  .476 

Permeability  <■  m  )  (various) 

Bulk  Thermal  Conductivity  o-f  Carbon  (  k^)  :  86.0  BtuZ-ft-hr — F 

Bulk  Specific  Heat  o-f  Carbon  (  Q,  )  :  .231  Btu/lbm-F 

Bulk  Density  o-f  Carbon  (Pc^^  70.3  lbm/-ft^ 

Thermal  Emissivity  o-f  Particles  (€  ):  .9 

Ambient  Temperature  (1^):  80.0  deg-F 

Ambient  Pressure  (  P ) :  14.7  psi 

Pressure  Di-f-f  erenti  al  (AP>  s  (various) 

Ambient  Oxygen  Concentration  (cfe,)  :  .0172  lbm/-ft 


82 


analysis  to  larger  values  of  Reynolds  number  would  exceed  the 
limitations  of  the  combustion  model.   Specifically,  the 
assumption  that  Darcy's  law  governs  the  internal  flow  for 

Reynolds  numbers  greater  than  5  would  be  questionable.   Simi- 

-9    2 
larly,  the  smallest  permeability  (m  =  .162  x  lo    f t  ) ,  at  a 

pressure  differential  of  .1  p.s.i.,  gave  an  average  Reynolds 
number  of  approximately  .005.   The  lower  limit  Reynolds  num- 
ber for  which  Darcy's  law  remains  valid  is  not  known.   There 
is  little  information  in  the  literature  that  addresses  a  lower 
limit  Reynolds  number.   Although  lower  Reynolds  numbers  may 
be  valid,  the  analysis  was  restricted  to  a  minimum  Reynolds 
number  of  .005.   As  can  be  seen  in  Figure  V.9,  the  "minimum 
initial  temperature"  for  sustaining  combustion  is  a  monotonically 
increasing  function  of  decreasing  slope  of  permeability.   In 
other  words,  as  permeability  increases,  the  temperature  needed 
for  sustained  combustion  also  increases.   Important  to  the  re- 
sults are  that  pore  velocity  also  increases  with  increasing 
permeability.   As  stated  previously,  insulated  boundaries 
were  imposed  on  the  porous  solid.   As  will  be  shown  in  Section 
V.E,  the  results  change  when  heat  transfer  occurs  at  the  sur- 
faces of  the  porous  solid. 
2.   Effects  of  Porosity 

Figure  V.IO  shows  the  dependence  of  the  "minimum 
initial  temperature"  for  sustaining  combustion  on  porosity. 
The  description  of  the  porous  medium  and  the  ambient  conditions 
used  in  the  following  analysis  is  shown  in  Table  V.4.   The 
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TABLE  V.4 Geometry  o-f  porous  medium  and  ambient  conditions 

for  the  porosity  analysis  in  Section  V/.D. 


Particle  Shape:  spherical 

Particle  Diameter  (d):  (various) 

Unit  Cell  Thickness  (Q  >!  (various) 

Spatial  Thickness  o-f  Porous  Medium  (L  ):       1.0  in 
Porosity  (  p  )  (various) 

Permeability  (m)  (various) 

Bulk  Thermal  Conductivity  o-F  Carbon  (  k^.)  :     86.0  Btu/-Ft-hr — F 

Bulk  Specific  Heat  o-f  Carbon  (  C^ )  :  .231  Btu/lbm-F 

Bulk  Density  of  Carbon  ( p^- >  :  70.3  Ibm/ft^ 

Thermal  Emissivity  of  Particles  (€):  .9 

Ambient  Temperature  (  X»)  :  SO.O  deg-F 

Ambient  Pressure  (  P ) :  14.7  psi 

Pressure  Differential  (AP) :  -.5  psi 

Ambient  Oxygen  Concentration  (<jfc>  )  :  .0172  Ibm/ft^ 
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results  show  that  as  porosity  increases,  the  temperature 
needed  for  sustaining  combustion  decreases  monotonically. 
This  behavior  is  opposite,  but  less  pronounced  than  that  ob- 
served for  the  permeability.   In  addition,  pore  velocity  de- 
creases as  porosity  increases.   It  seems  reasonable  that  an 
increase  in  porosity  provides  more  oxygen  per  volume  of  por- 
ous medium  for  combustion.   The  lowest  value  of  porosity  in 
the  analysis  (p  =  .48)  was  restricted  by  the  idealized  geometry 
shown  in  Figure  III.l.   For  the  idealized  geometry,  the  low- 
est porosity  occurs  for  the  ratio  of  particle  diameter  to 
unit  cell  thickness,  d/D,  equal  to  1.   The  maximiom  value  of 
porosity  in  the  analysis  was  limited  by  the  highest  porosity 
reported  for  spherical  particles  in  a  stable  configuration. 
As  stated  by  Scheidegger  [9],  this  value  is  .875. 
3.   Effects  of  Porous  Medium  Thickness 

Figure  V.ll  shows  the  dependence  of  the  "minimum 
initial  temperature"  for  sustaining  combustion  on  the  porous 
medium  thickness.   The  description  of  the  porous  medium  and 
ambient  conditions  used  in  the  analysis  is  shown  in  Table 
V.5.   The  range  of  thicknesses  (.25"-4.0")  investigated  was 
limited  (as  in  the  case  of  permeability)  by  the  Reynolds  num- 
ber.  As  seen  in  Figure  V.ll,  the  "minimum  initial  temperature" 
is  a  monotonically  decreasing  function  with  decreasing  nega- 
tive slope  of  porous  medium  thickness.   In  other  words,  as 
the  porous  medium  thickness  decreases,  the  temperature  needed 
to  sustain  combustion  increases.   In  addition,  the  pore  velocity 
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TABLE  V.5 Geometry  o-f  porous  medium  and  ambient  conditions 

for  the  thickness  analysis  in  Section  V.D. 


Particle  Shape:  spherical 

Particle  Diameter  ( d  ) :  (various) 

Unit  Cell  Thickness  (  D^s  (various) 

Spatial  Thickness  o-f  Porous  Medium  (  L  ):  (various) 

Porosity  (  p )  .476 

Permeability  (m)  (various) 

Bulk  Thermal  Conductivity  o-F  Carbon  (  k^)  :  86.0  Btu/ft-hr — F 

Bulk  Specific  Heat  of  Carbon  ( C^ )  :  .231  Btu/lbm-F 

Bulk  Density  of  Carbon  <0.):  70.3  Ibm/ft"^ 

Thermal  Emissivity  of  Particles  (€  ):  .9 

Ambient  Temperature  (  T^) :  80.0  deg-F 

Ambient  Pressure  (  F^) :  14.7  psi 

Pressure  Differential  (Ap> :  -.5  psi 

Ambient  Oxygen  Concentration  ^i>J '•  .0172  Ibm/ft 
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decreases  as  the  thickness  increases.   As  in  the  case  of  the 
permeability,  the  results  showed  a  wide  range  of  temperatures. 
The  "minimum  initial  temperature"  for  sustaining  combustion 
may  vary  as  much  as  700  degrees  Fahrenheit  over  a  range  of 
thicknesses  for  a  particular  porous  medium. 

E.   EFFECTS  OF  EXTERNAL  PARAMETERS 

A  similar  analysis,  as  discussed  in  the  previous  section, 
was  performed  for  the  pressure  differential.   The  boundary 
conditions  for  the  carbon  particle  temperature  equation  were 
then  changed  from  the  insulated  boundaries,  expressions  III.  4 5 
and  III. 46,  to  the  radiation  heat  transfer  boundary  conditions, 
expressions  III. 51  and  III.  52.   This  permitted  heat  transfer 
from  the  surfaces  of  the  porous  solid.   The  emissivity,  £,  was 
assiomed  equal  to  .9.   With  the  new  boundary  conditions,  the 
effects  of  permeability  and  pressure  differential  on  the  "mini- 
mum initial  temperature"  for  sustaining  combustion  were  again 
analyzed. 

1.   Effects  of  Pressure  Differential 

Figure  V.12  shows  the  effects  of  pressure  differential 
on  the  "minimum  initial  temperature"  for  sustaining  combustion. 
The  description  of  the  porous  medium  and  the  ambient  condi- 
tions used  in  the  following  analysis  is  shown  in  Table  V.6. 
The  range  of  pressure  differentials  considered  (.1-2.0  p.s.i.) 
was  governed  by  the  Reynolds  number  for  which  Darcy ' s  law  re- 
mains applicable.   This  point  was  discussed  in  Section  V.D. 
In  the  initial  work  of  Vatikiotis  [2] ,  the  maximum  pressure 
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TABLE  V.6 Geometry  o-f  porous  medium  and  ambient  conditions 

■for  the  pressure  di -f -Ferenti  al  analysis  in 
Section  V.E. 


Particle  Shape:  spherical 

Particle  Diameter  ( d  )  :  (various) 

Unit  Cell  Thickness  (Q^!  (various) 

Spatial  Thickness  o-f  Porous  Medium  (  L):  1-0  in 

Porosity  (  p  )  .476 

Permeability  (m)  (various) 

Bulk  Thermal  Conductivity  o-f  Carbon  (  k^)  :     S6.0  Btu/-ft-hr — F 

Bulk  Speci-fic  Heat  o-f  Carbon  (  C^  >  :              .231  Btu/lbm-F 

Bulk  Density  o-f  Carbon  (p^):  70.3  IbmZ-ft"^ 

Thermal  Emissivity  o-f  Particles  (€  ):  -9 

Ambient  Temperature  (T  ):  80.0  deg-F 

Ambient  Pressure  (  P ) :  14.7  psi 

Pressure  Di  f -ferenti  al  (AP)  :  (various) 

Ambient  Oxygen  Concentration  (Ob  )  :              .0172  Ibm/ft"^ 
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differential  considered  was  approximately  .35  p.s.i.   As  seen 
in  Figure  V.12,  the  "minimum  initial  temperature"  is  a  mono- 
tonically  increasing  function  of  pressure  differential.   The 
results  show  that,  depending  on  pressure  differential,  the 
minimum  temperature  for  sustaining  combustion  could  vary  as 
much  as  400  degrees  Fahrenheit  for  a  particular  porous  medium. 
In  addition,  pore  velocity  increases  as  pressure  differential 
increases. 

2.   Effects  of  Boundary  Conditions 

In  this  section,  the  results  of  allowing  heat  transfer 
at  the  boundaries  of  the  porous  solid  will  be  compared  to  the 
results  obtained  with  insulated  boundaries.   Figures  V.13  and 
V.14  show  the  results  of  the  permeability  and  the  pressure 
differential  analyses,  with  and  without  heat  transfer  from 
the  surface  of  the  porous  solid.   The  description  of  the  por- 
ous medium  and  ambient  conditions  for  the  permeability  and 
pressure  differential  analyses  are  shown  in  Tables  V.3  and  V.6, 
respectively.   As  stated  previously,  the  radiation  heat  trans- 
fer boundary  conditions  were  those  shown  by  expressions  III. 51 
and  III. 52.   As  seen  in  Figures  V.13  and  V.14,  heat  transfer 
from  the  boundaries  significantly  affects  the  "minimum  initial 
temperature"  at  the  lower  values  of  permeability  and  pressure 
differential  (or  at  lower  pore  velocities) .   The  results  show 
that  when  heat  transfer  occurs  from  the  boundaries,  the  "mini- 
mum initial  temperature"  is  no  longer  a  monotonic  function  of 
permeability  and  pressure  differential.   This  behavior  was 
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suggested  in  the  intial  work  of  Varikiotis  [2] ,  where  in  addi- 
tion to  radiation  heat  transfer,  convection  heat  transfer 
was  considered  at  the  boundaries.   The  additional  convection 
heat  transfer  causes  the  effects  produced  by  permitting  radia- 
tion heat  transfer  at  the  boundaries  to  become  more  pronounced. 
Moreover,  these  results  are  consistent  with  those  reported 
by  Vulis  [13] ,  Fontenot  [51] ,  and  Spalding  [52] .   The  reason 
for  this  behavior  was  summarized  in  Section  V.C  and  will  be 
explained  in  Section  V.G. 

3.   Effects  of  Initial  Conditions 

The  previous  analyses  of  determining  the  effects  of 
system  parameters  and  boundary  conditions  were  performed  using 
constant-value  initial  conditions  (i.e.,  slope  of  temperature 
distribution  was  equal  to  zero) .   The  constant-value  tempera- 
ture defined  the  "minimum  initial  temperature"  and  was  used 
as  a  basis  for  comparing  the  results.   To  also  show  the  effects 
of  the  initial  conditions  on  combustion,  a  single  porous  medium 
was  subjected  to  numerous  initial  conditions.   The  analysis 
was  performed  as  follows.   The  initial  conditions  of  the  car- 
bon and  air  temperatures  were  given  a  constant  slope  as  shown 
in  Figure  V.15.   With  the  slope  fixed,  the  temperature  distri- 
bution was  shifted  (i.e.,  moved  up  or  down)  until  the  thres- 
hold for  sustained  combustion  was  determined.   That  is,  shift- 
ing the  initial  temperature  distribution  slightly  lower  would 
result  in  the  porous  medium  cooling  to  ambient  temperature. 
This  procedure  was  repeated  for  several  initial  conditions  of 
temperature  having  different  constant  slopes  (both  negative 
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FIGURE  V.I5  Constant-slope  initial  condition. 


96 


and  positive) .   The  boundary  conditions  used  in  the  analysis 
were  the  Danckwerts  conditions,  expressions  III.  47  through 
III. 50,  for  the  air  temperature  and  oxygen  concentration,  and 
insulated  boundaries,  expressions  III. 4 5  and  III. 46,  for  the 
carbon  temperature.   The  results  of  the  analysis  are  shown  in 
Figure  V.16.   In  Figure  V.16,  the  average  temperature  is  plotted 
as  a  function  of  AT  defined  by  T(L)  -  T(0)  for  the  particular 
temperature  distribution.   The  average  temperature,  T    ,  is 
defined  as  (T(L)  +  T(0))/2.   As  seen  in  Figure  V.16,  the  curve 
is  not  symmetrical  about  AT  =  0 .   This  result  is  atrributed  to 
the  direction  of  the  pore  velocity.   The  largest  difference 
in  temperature  between  the  porous  solid  and  the  air  (and  great- 
est heat  transfer)  is  at  the  air  inlet  region  (i.e.,  x  =  0) . 
At  the  air  exit  region  (i.e.,  x  =  L),  the  temperature  differ- 
ence between  the  porous  solid  and  the  air  is  small.   The  com- 
bined effect  of  high  heat  transfer  at  the  air  inlet  region 
and  low  heat  transfer  at  the  air  exit  region  produces  the 
following  result.   There  is  a  greater  likelihood  that  sustained 
combustion  will  occur  for  a  positive  slope  temperature  dis- 
tribution than  for  a  negative  slope  temperature  distribution 
having  the  same  average  temperature  and  AT  as  defined  above. 
In  other  words,  for  a  positive  slope  initial  temperature  dis- 
tribution, sustained  combustion  occurs  at  a  higher  temperature. 

F.   EFFECTS  OF  PORE  VELOCITY 

The  way  the  parameters  (i.e.,  permeability,  porosity, 
pressure  differential,  and  porous  medium  thickness)  affected 
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combustion  is  in  proportion  to  their  effect  on  pore  velocity. 
In  other  words,  the  results  of  the  previous  sections  showed 
that  pore  velocity  was  affected  by  the  parameters  in  the  same 
manner  as  the  "minimum  initial  temperature"  for  sustained  com- 
bustion (i.e.,  minimum  temperature  increased  as  pore  velocity 
increased,  and  decreased  as  pore  velocity  decreased) .   An 
exception  to  this  occurred  when  heat  transfer  was  permitted 
from  the  boundaries.   This  suggests  that  pore  velocity  is  the 
dominant  variable  affecting  the  combustion  process.   Since 
pore  velocity  is  a  function  of  the  parameters  analyzed  (i.e., 
permeability,  porosity,  pressure  differential,  porous  medium 
thickness) ,  it  seems  reasonable  that  the  "minimum  initial 
temperature"  for  sustained  combustion  may  be  posed  as  a  func- 
tion of  pore  velocity.   With  this  approach,  the  "minimum  initial 
temperature"  for  sustained  combustion  of  a  porous  medium  with 
unknown  properties  could  be  determined  by  specifying  the 
magnitude  of  the  pore  velocity.   Moreover,  pore  velocity  can 
be  represented  by  specific  volumetric  flowrate  of  filter  velocity 
(i.e.,  V/A)  which  is  easily  measured.   Figure  V.17  shows  the 
"minimum  initial  temperature",  defined  in  Section  V.D,  as  a 
function  of  the  filter  velocity  at  ambient  temperature.   The 
graph  was  made  from  the  data  obtained  in  the  analyses  of  Sec- 
tions V.D  and  V.E  (for  insulated  boundaries  on  the  porous 
solid) .   The  data  in  Figure  V.18  has  a  spread  of  approximately 
±100  degrees  Fahrenheit.   The  equation  for  the  curve  shown 
in  Figure  V.17  is. 
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Tg^   =   100  ln^(Q/10)  +  650  (V.l) 


where  T   represents  the  "minimum  initial  temperature"  dis- 
tribution  in  degrees  Fahrenheit,  and  Q  is  the  filter  velocity 
in  ft/hr  at  an  ambient  temperature  of  80  degrees  Fahrenheit. 
Figure  V.17  shows  promise  as  a  simple  method  for  determining 
the  "minimum  initial  temperature"  for  which  sustained  combus- 
tion can  be  expected  in  a  porous  medium.   Perhaps  by  refining 
the  formulation  for  calculating  the  pore  velocity,  there  would 
be  a  smaller  spread  in  the  data. 

G.   ANALYSIS  OF  RESULTS 

The  Semenov  model,  discussed  in  Section  III.C  and  illus- 
trated in  Figure  III. 2,  may  be  used  to  explain  the  results 
reported  in  Sections  V.B  through  V.F.   Before  continuing,  cer- 
tain aspects  of  the  Semenov  model  must  be  emphasized.   First, 
the  Semenov  model  does  not  govern  combustion.   It  is  an  analy- 
tical tool  which  is  used  to  describe  the  behavior  of  combustion 
Secondly,  the  Semenov  model  as  represented  in  Figure  III. 2 
is  for  a  point  in  the  porous  medium.   A  "Semenov  surface" 
would  better  explain  combustion  in  the  porous  medium.   This 
will  be  illustrated  below.   It  follows  from  the  Semenov  theory, 
as  discussed  by  Vulis  [13],  that  "ignition"  conditions  (i.e., 
sustained  combustion)  are  determined  by  all  the  conditions  of 
the  combustion  problem  in  a  system.   This  is  to  say,  that  any 
change  in  system  parameters  (e.g.,  pressure  differential, 
initial  or  boundary  conditions)  for  a  given  porous  medium 
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would  result  in  a  different  Semenov  surface.   Lastly,  the 

Semenov  model  is  based  on  quasi-steady  results.   For  example, 

the  convection  heat  transfer,  q„ ,  and  the  heat  generation, 

R  ,  will  be  equal  at  their  intersection  (see  Figure  III. 2). 

Conversely  for  the  combustion  model,  at  a  particular  point 

and  time,  q,  may  be  smaller  or  larger  than  R  depending  on 
X,  g 

whether  the  temperatures  are  increasing  or  decreasing.   How- 
ever, if  the  rates  of  change  in  the  temperatures  are  small, 
the  Semenov  model  can  be  used  to  analyze  the  results  of  the 
combustion  model . 

Figure  V.18  illustrates  a  "Semenov  surface"  generated 
from  the  results  of  the  last  example  in  Section  V.C  (Figures 
V.7  and  V.8) .   In  Figure  V.18,  R   is  a  heat  generation  sur- 
face, and  represents  the  particle  temperature-heat  generation 
relationship  during  the  problem.   The  q.  surface  represents 
the  convection  heat  transfer  as  a  function  of  particle  and 
air  temperature  (expression  III. 14).   The  particular  heat 
transfer  surface  in  Figure  V.i8  represents  the  instant  the 
point  at  x/L  =  0  transitioned  from  the  kinetic  regime  to 
the  diffusion  regime  of  combustion.   This  can  be  seen  by  the 
tangent  point,  I,  between  the  heat  generation  surface  and  the 
heat  transfer  surface.   A  difference  between  "classical" 
Semenov  S-curves  (Figure  III. 2)  and  those  shown  in  the  heat 
generation  surface  of  Figure  V.18  is  the  sudden  drop  in  heat 
generation  once  a  certain  temperature  has  been  reached.   This 
is  only  observed  in  the  region  away  from  the  air  inlet  (x/L  =  0) 
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FIGURE  V.18   Semenov  surface. 
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and  is  attributed  to  the  total  depletion  of  oxygen  in  that 
region  of  the  porous  medium.   One  might  say  this  is  a  weak- 
ness in  the  Semenov  model  when  applied  to  porous  media. 
Conversely,  it  might  be  argued  that  when  the  oxygen  is  totally 
depleted  the  combustion  problem  becomes  a  heat  transfer  prob- 
lem (in  a  pure  sense)  of  which  the  Semenov  model  has  no 
application.   Also  in  the  region  where  the  oxygen  has  been 
depleted,  it  is  not  necessary  for  q^  to  intersect  R  . 

In  order  to  visualize  the  results  of  the  combustion  model 
more  clearly,  the  "Semenov  surface"  is  abandoned  for  Semenov 
curves  (illustrated  in  Figure  III. 2).   Figures  V.19,  V.20, 
and  V.21  show  that  Semenov  model  representations  for  locations 
x/L  =  0,  x/L  =  .5,  and  x/L  =  1,  respectively,  of  example  C  in 
Section  V.C  (transient  results  shown  in  Figures  V.7  and  V.8). 
The  heat  generation  curves  are  shown  along  with  heat  transfer 
lines  representing  the  film  cooling  at  two  different  points 
in  time.   During  this  discussion,  the  term  "film  cooling"  will 
be  used  to  mean  the  heat  transfer  from  the  particle  to  the 
air.   "Convection  heat  transfer"  will  refer  to  the  transport 
of  energy  by  the  internal  flow.   In  addition  to  film  cooling, 
heat  transfer  by  conduction  and  radiation  occurs  within  the 
porous  medium.   This  has  an  effect  of  slightly  curving  the 
heat  transfer  line  upward  (i.e.,  heat  transfer  becomes  a  non- 
linear vice  a  linear  function  of  temperature) .   However  in 
this  problem,  calculations  showed  that  the  ratio  of  film 
cooling  to  conduction  and  radiation  was  approximately  5 ,  and 
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therefore,  conduction  and  radiation  will  not  be  considered 
in  order  to  simplify  the  graphical  analysis. 

Comparing  Figures  V.19-V.21  shows  that  the  maximum  heat 
generation  rates  are  different  depending  on  the  location  in 
the  porous  medium.   It  was  mentioned  in  the  general  discussion 
of  the  Semenov  model  (Section  III.C)  that  implicit  in  the 
shapes  of  the  S-curves  are  the  effects  of  oxygen  concentra- 
tion.  In  other  words ,  decreasing  the  oxygen  supply  will 
cause  the  S-curves  to  peak  at  lower  temperatures.   Since  it  is 
more  difficult  for  the  oxygen  to  penetrate  as  the  depth  of 
the  porous  medium  increase  (i.e.,  oxygen  is  consumed  as  the 
flow  passes  through  the  porous  medium) ,  it  is  a  reasonable 
outcome  that  the  maximum  heat  generation  rates  decrease  as 
the  depth  increases.   The  difficulty  for  the  oxygen  to  pene- 
trate as  the  depth  increases  produces  another  effect  as  seen 
in  Figure  V.21.   Unlike  those  in  Figures  V.19  and  V.20,  the 
S-curve  and  heat  transfer  curve  does  not  have  a  tangent  point 
which  defines  the  "critical  ignition  condition".   This  be- 
havior is  called  noncritical  combustion,  and  as  discussed  by 
Vulis  [13] ,  is  characteristic  of  oxygen  poor  combustion. 
Moreover,  these  observations  are  supported  by  the  experimental 
results  of  Thomas,  Stevenson,  and  Evans  [15].   In  their 
experiments,  it  was  shown  by  decreasing  the  ambient  partial 
pressure  of  oxygen,  that  a  critical  combustion  process  (i.e., 
one  in  which  a  "temperature  jump"  occurs)  transforms  to  a 
noncritical  combustion  process. 
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As  seen  in  Figure  V.7,  the  temperature  first  rises  in  a 
local  region  of  the  porous  medium  (in  this  case  near  x/L  =  1) 
before  moving  as  a  wave  toward  the  air  inlet  surface.   The 
direction  of  movement  was  a  characteristic  of  all  the  prob- 
lems in  which  sustained  combustion  was  observed.   This  behavior 
is  explained  as  follows.   Temperature  rises  in  a  region  where 
the  heat  generation  is  greater  than  the  heat  transfer  between 
the  particles  and  the  air.   The  region  of  higher  temperature 
widens  as  heat  is  transferred  by  conduction  and  radiation. 
In  addition,  the  increasing  temperature  simultaneously  in- 
creases the  reaction  rate.   At  some  moment,  the  oxygen  enter- 
ing through  the  upstream  side  of  the  high  temperature  region 
is  totally  consumed  within  that  region.   When  this  occurs, 
the  simultaneous  increase  in  temperature  and  reaction  rate 
will  occur  only  on  the  upstream  side  of  the  region  of  high 
temperature.   As  the  increasing  temperature  moves  towards  the 
air  inlet  surface,  the  front  of  the  region  of  depleted  oxygen 
also  moves  nearer  the  air  inlet  surface.   This  results  in  the 
temperature  response  appearing  to  move  in  the  form  of  a  wave 
toward  the  air  inlet  surface.   It  is  interesting  to  note,  that 
the  magnitude  of  heat  transfer  by  conduction  and  radiation  in 
the  increasing  temperature  region  was  approximately  twice 
that  of  the  convection  heat  transfer  at  the  outset  of  the 
example  (Figures  V.7  and  V.8).   It  seems  reasonable  that  by 
increasing  the  pore  velocity  to  a  sufficient  magnitude,  the 
region  of  increasing  temperature  would  not  propagate  but  blow 
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out  (i.e.,  toward  the  x/L  =  1  surface) .   In  fact,  this  is 
observed  for  the  problems  in  which  the  porous  medium  cools 
to  ambient  temperature. 

Up  to  this  point,  the  thermal  response  of  the  porous  medium 
has  been  discussed  in  terms  of  "sustained  combustion"  or 
"cooling  to  ambient  temperature".   The  Semenov  curves  gener- 
ated from  the  results  of  the  combustion  model  can  be  used  to 
determine  transition  from  kinetic  to  diffusion  combustion. 
The  kinetic  and  diffusion  regimes  of  combustion  are  discussed 
in  Section  III.C.   In  a  practical  sense,  the  exothermic  proc- 
ess occurring  in  the  kinetic  regime  is  synonymous  with  the 
term  "oxidation".   The  process  occurring  in  the  diffusion 
regime  is  simply  called  "combustion".   From  everyday  experi- 
ence, one  knows  that  oxidation  can  be  a  slow  process,  occurring 
at  low  reaction  rates.   Conversely,  combustion  can  produce 
considerable  heat,  reflected  by  high  temperatures  and  reaction 
rates.   As  will  be  shown,  transition  from  oxidation  to  combus- 
tion occurs  rapidly.   In  addition,  the  porous  medium  may  be 
undergoing  oxidation  in  one  part,  and  combustion  in  the  other. 
It  is  important  to  understand  the  behavior  of  the  thermal 
process  at  this  level.   Especially,  if  there  is  a  need  to 
closely  control  the  process. 

Each  point  within  the  porous  medium  will  transition  from 
oxidation  to  combustion  at  different  times.   This  can  be 
seen  in  Figures  V.19-V.21  which  were  obtained  from  the  results 
of  example  C.   The  locations  x/L  =  1,  x/L  =  .5,  and  x/L  =  0 
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transitioned  at  12.5  seconds,  22  seconds,  and  41  seconds, 
respectively.   Therefore  at  the  start  of  the  problem,  the 
porous  medium  as  a  whole,  was  undergoing  oxidation.   More- 
over, transition  from  oxidation  to  combustion  started  at 
x/L  =  1  and  moved  to  x/L  =  0  with  time.   This  is  consistent 
with  the  behavior  observed  for  the  temperature  response.   As 
discussed  by  Vulis  [13] ,  the  theoretical  value  for  "critical 
ignition"  temperature  of  the  carbon  can  be  calculated  using 
the  "N,  N.  Semenov  equation"  given  by. 


^I   =  W~^^   '    (1- 4R^T^/E)^/^]  (V.2) 


where  E  is  the  activation  energy  of  the  Arrhenius  expression, 

R   is  the  universal  gas  constant,  and  T  is  the  absolute  air 
u  ^  a 

temperature.   This  equation  is  derived  by  equating  the 
Arrhenius  expression  III. 16  to  the  convection  heat  transfer 
expression  III. 14,  and  requiring  the  slopes  be  equal.   Note 
that  the  theoretical  "critical  ignition"  temperature  is  a 
function  of  air  temperature  and  activation  energy  only.   Since 
the  transition  from  oxidation  to  combustion  at  x/L  =  1  shows 
a  noncritical  behavior,  equation  V.l  does  not  apply.   The 
theoretical  "critical  ignition"  temperatures  for  locations 
x/L  =  0  and  x/L  =  .5  are  1831  degrees  Fahrenheit  and  1492 
degrees  Fahrenheit,  respectively.   The  values  obtained  by  the 
combustion  model  (Figures  V.19  and  V.20)  are  approximately 
1825  degrees  Fahrenheit  for  x/L  =  0,  and  1440  degrees  Fahren- 
heit for  x/L  =  .5.   Once  the  "critical  ignition"  temperature 
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was  reached,  the  increase  to  the  "combustion"  temperature 
occurred  in  less  than  1  second  for  both  cases .   Examples  A 
and  B  in  Section  V.C  (Figures  V.2-V.6)  were  the  same  problem  as 
example  C,  discussed  here,  except  combustion  was  initiated 
differently.   Example  C  was  initiated  with  the  porous  medium 
at  a  constant  temperature  known  to  sustain  combustion.   In 
examples  A  and  B,  the  porous  medium  was  subjected  to  a  heat 
flux  at  the  air  inlet  surface  (x/L  =  0)  with  the  porous  medium 
at  ambient  temperature.   The  amount  of  time  for  which  the  heat 
flux  was  applied  determined  if  sustained  combustion  occurred. 
Sustained  combustion  was  observed  for  example  A  in  Section 
V.C.   The  amount  of  time  the  heat  flux  was  applied  in  example 
B  was  not  sufficient  for  sustaining  combustion  and  the  porous 
medium  cooled  to  ambient  temperature.   It  is  interesting  to 
note  that  the  temperature  at  x/L  =  0  in  example  A  reached  the 
"critical  ignition"  temperature  of  1825  degrees  Fahrenheit. 
The  maximum  temperature  obtained  at  x/L  =  0  in  example  B  was 
approximately  1650  degrees  Fahrenheit.   There  is  good  agreement 
between  the  theoretical  values  of  "critical  ignition"  tempera- 
ture and  those  obtained  by  the  combustion  model . 

The  pore  velocity  significantly  affects  the  combustion 
process.   This  was  discussed  in  Section  V.F.   Increasing  pore 
velocity  has  two  effects.   First,  the  internal  heat  transfer 
coefficient,  h,  is  increased,  thereby,  increasing  the  amount 
of  film  cooling.   Secondly,  the  supply  of  oxygen  is  made  greater 
(i.e.,  convection  mass  transfer  of  oxygen  molecules  increases). 
The  combined  effect  is  that  by  increasing  pore  velocity,  higher 
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local  temperatures  are  needed  for  transition  of  the  porous 
solid  from  oxidation  to  combustion.   In  other  words,  assume 
that  transition  occurs  at  a  particular  air  temperature  for  a 
given  pore  velocity.   Increasing  the  pore  velocity  without  an 
increase  in  air  temperature  may  not  be  sufficient  for  transi- 
tion to  occur.   This  is  illustrated  by  Figure  V.22.   Increas- 
ing the  pore  velocity  will  increase  the  slope  of  the  heat 
transfer  line,  and  also  change  the  shape  of  the  heat  genera- 
tion curve  as  shown.   It  is  apparent  that  an  increase  in  air 
temperature  is  needed  to  make  the  steeper-sloped  heat  transfer 
line  tangent  to  the  heat  generation  curve.   In  addition,  higher 
pore  velocities  produce  higher  combustion  temperatures  once 
transition  occurs.   In  the  general  discussion  of  the  Semenov 
model,  it  was  pointed  out  that  the  temperature  dominates  the 
combustion  process  in  the  kinetic  regime.   Increasing  the 
oxygen  supply  will  have  the  greatest  effect  in  the  diffusion 
regime  (i.e.,  increasing  the  combustion  temperature).   This 
behavior  described  by  the  Semenov  model,  combined  with  the 
effects  produced  by  convection  heat  transfer  (energy  trans- 
port by  internal  flow)  and  to  a  lesser  extent,  heat  transfer 
by  conduction  and  radiation,  will  determine  if  the  porous 
medium  will  undergo  sustained  combustion  for  an  increase  in 
pore  velocity. 

Figure  V.23  shows  the  effects  of  decreasing  the  pore 
velocity.   The  slope  of  the  heat  transfer  line  becomes  smaller, 
and  the  heat  generation  curve  changes  shape  as  shown.   In 
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particular,  lower  combustion  temperatures  will  occur  in  the 
diffusion  regime.   Unlike  the  increasing  pore  velocity  case 
above,  there  is  less  certainty  in  predicting  the  behavior  of 
combustion  when  pore  velocity  is  decreased.   It  appears  that 
conduction  and  radiation  heat  transfer  become  more  of  a  factor 
at  the  lower  pore  velocities.   This  seems  reasonable  since  the 
ratio  of  conduction  and  radiation  heat  transfer  to  convection 
heat  transfer  increases  with  decreasing  pore  velocity.   In 
addition,  the  effects  of  heat  transfer  from  the  boundaries  at 
the  lower  pore  velocities  must  also  be  considered  (shown  in 
Figure  V.13  and  V.14).   The  combined  effects  of  all  the  heat 
transfer  mechanisms,  including  the  boundary  conditions,  will 
determine  if  the  porous  medium  will  undergo  sustained  combus- 
tion as  the  pore  velocity  is  decreased. 

Lastly,  the  effects  of  the  boundary  conditions  must  be 
considered  separately.   The  results  of  changing  the  boundary 
conditions  were  presented  in  Section  V.E.   It  was  observed 
that  for  heat  transfer  occurring  at  the  surfaces  of  the  porous 
solid  Ce.g.,  radiation  boundary  conditions),  higher  initial 
temperatures  were  needed  to  sustain  combustion.   Allowing  heat 
transfer  from  the  surfaces  of  the  porous  solid  (either  by 
convection  or  by  radiation  or  by  both) ,  affects  the  combustion 
process  in  the  boundary  condition  regions  as  follows.   Return- 
ing to  the  Semenov  model,  there  is  nothing  explicitly  asso- 
ciated with  the  boundary  conditions  that  would  change  the  shape 
of  the  heat  generation  curve.   In  addition,  since  pore  velocity 
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determines  the  magnitude  of  the  internal  heat  transfer  coeffi- 
cient, it  appears  that  the  slope  of  the  heat  transfer  line 
will  also  be  unaffected.   However  up  to  this  point,  only  film 
cooling  was  considered  in  the  heat  transfer  line  of  the  Semenov 
model.   Noting  that  heat  transfer  from  the  porous  solid  boun- 
daries results  in  steeper  temperature  gradients  in  regions 
near  the  boundaries,  heat  transfer  by  conduction  and  radiation 
must  also  be  considered.   With  the  additional  components  of 
heat  transfer,  the  heat  transfer  line  will  curve  upward.   This 
means  that  the  heat  transfer  changes  from  a  linear  function 
to  a  nonlinear  function  of  temperature.   The  change  of  the 
shape  is  illustrated  in  Figure  V.24.   As  can  be  seen,  the 
increased  heat  transfer  by  conduction  and  radiation  may  pre- 
vent transition  from  oxidation  to  combustion.   This,  in  turn, 
will  restrict  the  heat  generation  to  the  kinetic  regime.   It 
is  reasonable  to  expect  that  as  the  reaction  moves  towards  the 
air  inlet  boundary,  the  combined  effects  of  greater  heat  trans- 
fer and  lower  heat  generation  will  oppose  sustained  combustion. 
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VI.   CONCLUSIONS 

The  main  objective  of  developing  a  model  which  predicts 
the  combustion  behavior  of  a  porous  medium  has  been  achieved. 
The  results  of  the  mathematical  model  are  in  good  agreement 
with  those  obtained  by  experimental  methods.   Moreover,  the 
analyses  presented  have  shown  that  combustion  and  heat  trans- 
fer in  porous  media  is  a  complex  process  involving  the  inter- 
action of  heat  transfer  and  heat  generation.   This  interaction 
depends  upon  the  geometric  parameters  of  a  porous  medium, 
boundary  conditions,  environmental  and  initial  conditions. 
The  behavior  of  a  system  can  change  radically  by  altering  any 
number  of  parameters.   This  point  was  demonstrated  by  examples 
A  and  B  in  Section  V.C  where  a  difference  of  one  second  in 
applying  a  surface  heat  flux  meant  combustion  or  extinguish- 
ment.  The  results  have  also  shown  that  conduction  and  radia- 
tion between  particles  may  play  a  significant  role  in  deter- 
mining combustion  behavior  and  should  be  accounted  for.   The 
computer  program  makes  it  possible  to  look  at  a  large  number 
of  cases.   Similar  analyses  by  experimental  methods  would  be 
economically  impractical. 

Based  on  literature  surveys   performed  during  this  inves- 
tigation, it  appears  that  much  of  the  engineering  develop- 
ment associated  with  porous  media  invoJ-ves  a  trial  and  error 
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process.   It  seems  reasonable  that  an  analytical  model  would 
be  an  essential  tool  for  the  engineer  to  employ  in  the  de- 
sign process.   In  this  investigation,  the  model  has  been 
assessed  for  its  applicability,  and  to  some  extent,  its 
accuracy.   It  is  hoped  that  in  the  future  the  model  will  be 
used  to  determine  the  effects  of  the  design  variables  (i.e., 
the  geometric  parameters,  pressure  differential,  etc.)  on 
performance.   Specifically,  the  model  shows  promise  for 
evaluating  the  combustion  efficiency  and  stability  of  a  sys- 
tem.  It  is  in  this  capacity  that  the  combustion  and  heat 
transfer  model  will  be  of  greatest  value. 
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APPENDIX  A 
FORMULATION  OF  FIELD  EQUATIONS 

1.   PRESSURE  DISTRIBUTION  EQUATION 

Darcy • s  law  for  one  dimensional  flow,  neglecting  body 
forces,  is. 


Q   =   -  -(|^)  (A.l) 

U  dx 


Substituting  in  the  Dupuit-Forcheimer  relation,  and  solving 
for  u,  equation  A.l  becomes. 


The  continuity  equation  (derived  in  Appendix  B)  is. 


3(pp^)    . 


Substituting  equation  A. 2  into  equation  A. 3  yields. 


3(pp  )       mp   , 


Expanding  terms,  equation  A. 4  becomes. 


4.  <i^.i|S.l  3P,|P.^!4!£i=   0       (A.5) 
dx      '^a       m  3x    y  8x  dx    mp     3t 
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Equation  A. 5  with  associated  boundary  conditions  (presented  in 
Section  II. B)  is  solved  for  the  pressure  and  pressure  gradient 
distribution.   The  pressure  gradient  is  then  substituted  into 
equation  A. 2  to  obtain  the  pore  velocity. 

2.   POROUS  SOLID  HEAT  TRANSFER  EQUATION 

To  perform  energy  balances  on  both  the  porous  solid  and 
on  the  air,  a  differential  volume  of  porous  medium  was  segre- 
gated into  respective  volumes  of  constituents,  that  is, 
dV  =  (l-p)dV  for  the  solid,  and  dV  =  pdV  for  the  air  (shown 
in  Figure  A.l).   The  convention  used  for  the  energy  balance 
of  an  arbitrary  differential  volume,  dV,  is. 


Heat  into 
dV 


Heat 
Generation 


Heat  out 
of  dV 


Increase  in 
Internal  Energy 


The  heat  transfer  mechanisms  considered  for  the  carbon 
particles  are  conduction,  radiation  heat  transfer  between 
particles,  convection  heat  transfer  from  the  particles  to 
the  air,  and  heat  generation.   Applying  the  above  convention, 
the  energy  balance  on  a  differential  volume  of  porous  solid 
is. 


^l-P^^cond^^l   ^  (l-P)q^ad^A|   +  q    dA' 
X  X     ^ 


(A. 6) 


(1-P^^cond^^ 


x+dx 


-K  (l-p)q^^^dA 


^ ,   +  q     dA' 
x+dx    ^conv 


*  (i-pjqint-^v 
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FIGURE  A.l  Separating  a  differential  volume  of 

porous  medium  into  respective  volumes 
of  solid  and  air. 
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Representing  terms  on  the  right  side  of  expression  A. 6  by 
Taylor  series  expansions  (neglecting  higher  order  terms) , 
the  energy  balance  becomes , 


(l-p)q^O„d'^A    +  (l-p)q^^^dA    +  q    dA- 


X 


X 


(A. 7) 


=   ^^"P^^ond^   -^  37^^1-P^^cond^^^^  ^  (l-P)q..^dA 

X 


^rad' 


X 


+  ^[(l-p)q   ,]dxdA  +  q     dA'  +  (l-p)q.  .dV 
3x    ^  ^rad        ^conv         ^  ^int 


Subtracting  terms,  and  rearranging,  expression  A. 7  becomes. 


-  I— [(l-p)q    ,]dV  -  -|-[(l-p)q   ,]dV  -  q     dA' 
dx  ^  ^cond'      ax"-    ^  ^rad      ^conv 


+  q    dA'   =   (l-p)q.  ^dV 
^gen  ^'^mt 


(A. 8) 


Substituting  the  following  expressions  into  equation  A. 8, 


dT 


^cond 


-  k 


e  3x 


Fourier's  law 


(A. 9) 


dT 


^rad 


=   -  k 


r  9x 


Radiation  analogy  to      (A. 10) 
Fourier's  law 


^conv  ^   ^^"^c  "  ^a^      Newton's  law 


(A. 11) 


q      =   R 
^gen      g 


Heat  generation 


(A. 12) 
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^int  ^  ^c  ^c  TT       Internal  energy  (A.  13) 


yields , 


|-[(l-p)  (k  +k  )  -5-^]dV  -  h(T  -T^)dA'  +  R  dA'  (A. 14) 

3T 


Dividing  through  by  dV,  and  defining  dA'/dV  as  z,    the  specific 
internal  area  (i.e.,  surface  area  per  unit  volume),  equation 
A.  14  becomes. 


ijjKl-pXk^+k^)  if]  -  hz(T^-T^)  +  R^z  (A.15) 

3T 
=   (1-p)  P^  c^  ^ 


The  expressions  used  to  obtain  the  values  of  the  properties 
and  parameters  in  equation  A.15  are  presented  in  Section 
III.E. 

3.   AIR  HEAT  TRANSFER  EQUATION 

The  formulation  of  the  air  heat  transfer  equation  will 
begin  with  the  general  one  dimensional  energy  equation. 


9T 
PPait^^^l-')   =   fe^P^alir)  ^  ^^^V^a^  ^  ^x       ^^'^'^ 


l-(puP)  -|-(p  U  T    ) 

3x  ^  dx   ^        XX 
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As  in  the  porous  solid  heat  transfer  equation,  the  time  and 
position  dependent  porosity  appears  inside  the  differential. 
Expanding  terms  and  neglecting  body  forces,  the  air  energy 
equation  becomes , 


-   T— (puP)  -  TT— (put   ) 
9x  ^        8x  '^    XX 


Consider  the  momentum  equation  for  the  x-direction  (neglecting 
body  forces) , 


It^PPa^^   =   -  fe^PPa^^)  -  Is^^P-^xx^  -  fe^PP)         (^-^^^ 


and   the  continuity  equation. 


It^PPa^      =     -^k^PPa^    "   PPal^  ^^'^^^ 


Multiplying  the  continuity  equation  through  by  u,  and  substi- 
tuting this  into  equation  A. 18,  the  momentum  equation  becomes. 


P^al^  =  -  PPa-0-  k'P^'  -k<P^xx>         '^-2°' 


Multiplying  equation  A. 20  by  u  and  noting  that. 


^u  8u  ^       2  3u  ,^    ^,  , 

upp^  Dt   =   ^PPa  3t  -^  PPa^  T^  ^^-21) 
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equation  A. 20  becomes, 


^P  Pa  ST  =   -  ^  ll^^PP^  -  ^  k^P^xx^  ^^-22) 


The  energy  equation  A. 17,  after  substituting 

1  2 
D(yu  ) /Dt  =  uDu/Dt  and  expanding  terms,  is, 

ST 


„  8u      8  pP 
P  P  7T u    '^ 

^   8x       3x 


3  -     ,         3u 

U  TS — (PT^    )   -  P  T     :r— 

3x  '^   XX     ^  XX  3x 


Substituting  equation  A. 22  into  the  above  energy  equation  and 
cancelling  terms,  equation  A. 23  becomes. 


3T 


The  viscous  dissipation  term,  px  3u/3x,  is  neglected  because 
the  fluid  is  a  gas  flowing  at  a  low  velocity.  Therefore,  the 
energy  equation  for  the  air  in  the  porous  medium  is, 

PP  ^  =   |-(pl^  -^)    +  hz(T  -T  )  -  pp|^        (A. 25) 
^^aDt     dx  '^    a    dx  ca    ^3x 

With  specific  enthalpy  for  a  gas  defined  by, 

k      =  e  +  P/p  (A. 26) 
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De/Dt  can  be  expressed  as. 


De   _   Dfi    1  DP  .  P  ^a  n,    ni\ 

Dt   "   Dt  "  p^  Dt     2  Dt  KP^.^n 

^  Pa 

2 

Multiplying  the  continuity  equation  through  by  P/(pp  )  gives, 

a 


p 

°Pa 

P      DP 

P     3u 

Pa 

Dt 

PPa    ^^ 

p^    3x 

(A. 28) 


Substitution  of  the  above  expression  into  expression  A. 27, 
the  substantial  derivative  of  internal  energy  can  be  expressed 
as. 


De      ^      D^        J^  DP        ^  iii        _l_  DP  ,,    29) 

Dt  Dt   ~   p^    Dt   "    p^    9x   "    pp^   Dt  ^A.^^; 

a  a  a 


Substituting  expression  A. 29  into  equation  A. 25,  the  energy 
equation  reduces  to. 


PPa5|-P§|-^^   =   k'P'^a^'  ^'^^'(V^a'       ^^-^"^ 


The  specific  enthalpy  can  be  represented  in  terms  of  tempera- 
ture ,  as , 


dfi   =   T  ds  +  -  dP  (A.  31) 


where  s  is  specific  entropy.   For  a  perfect  gas,  ds  may  be 
expressed  by. 
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■^^  =  =p  ¥-  ^<5P  l^-^^' 


Substituting  expression  A. 32  into  expression  A. 31,  and  can- 
celling terms,  enthalpy  can  now  be  written  as. 


dk      =   c   dT  (A.  33) 

p 


or. 


r-.t  DT 

^     =      c      -^  (A. 34) 

Dt       a  Dt 


Substituting  expression  A. 34  into  equation  A. 30,  the  air 
energy  equation  or  heat  transfer  equation  becomes. 


3  T  9T 


Initial  results  showed  that  pP  changed  slowly  with  time. 
Therefore,  the  substantial  derivative  of  pP ,  as  shown  here. 


§t(pp)    =    It^p^^  ■"  ^  Is^^p^^  ^^'^^^ 


is  reduced  to  u8(pP)/8x.   The  expressions  used  to  obtain  the 
properties  and  parameters  in  the  coefficients  of  equation 
A. 35  are  presented  in  Section  III.E. 
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4.   OXYGEN  MOLECULE  DIFFUSION  EQUATION 

The  final  consideration  in  formulating  the  field  equa- 
tions for  the  model  is  the  transport  of  oxygen  molecules. 
The  oxygen  molecule  transport  equation  is  obtained  by  a  con- 
servation of  species  belance  on  the  differential  volume  of 

air,  dV  =  pdV.   The  convention  used  for  the  species  balance 
a 

is  given  by. 


0^  into 
^   dV 


O2  out  ^02        +     O2 

of  dV    Consumption   Accumulation 


The  transport  mechanisms  considered  were  diffusion  resulting 
from  concentration  gradients,  convection,  and  consumption  of 
oxygen  by  combustion.   Applying  the  above  convention,  the 
species  balance  on  the  oxygen  becomes. 


P  nij  .  £j-  dA   +  p  m     dA    =   p  m ,  .  -  _  dA 
^    diff        '^  conv  ^    diff 

X  X 


(A. 37) 


x+dx 


+  p  m     dA 
^    conv 


x+dx 


+  m     dA '  +  p  m    dV 
cons       ^  ace 


Representing  the  terms  on  the  right  side  by  Taylor  series 
expansions  (neglecting  higher  order  terms)  ,  the  species  balance 
becomes , 


pin,  .^^dA   +  pm     dA    =   p  m ,  .  ^^  dA 
^  diff        ^  conv  ^  diff 

X  X 


(A. 38) 


+  T:r(pm, .  ^^)dxdA  +  pm    dA   +  ^r— (pm     )  dxdA  +  m    dA'  +pm   dV 
3x  '^  diff        ^   conv      3x  ^  conv'         cons     ^  ace 

X 
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Cancelling  terms  and  rearranging,  equation  A. 38  becomes, 


k<P™diff'<^^  -  b'P"conv'^^  -  "cons"^^'   =  P-^acc"^^   '^-^S' 


Substituting  the  following  expressions  into  equation  A. 39, 


diff       e  3x 


Pick's  law 


(A. 40) 


^conv  =   ^  ^ 


Convection  transport      (A. 41) 


"'cons   =   ^0. 


Consumption  by 
combustion 


(A. 42) 


m 


ace 


at 


Accumulation 


(A. 43] 


yields. 


k'P  "elt'^^  -  kc^P*'^^  -  \'^-      =   P  H  <^^        <^-44) 


Dividing  both  sides  by  dV,  and  letting  dA'/dV  equal  the 
specific  internal  area,  z,  the  oxygen  molecule  diffusion 
equation  becomes , 


8_ 
8x 


(P-'eH)  -|3r<-P*'  -  ^0^^   =   P|! 


(A. 45) 


The  methods  and  expressions  for  obtaining  the  properties  and 
parameters  in  the  coefficients  of  equation  A. 45  are  presented 
in  Section  II. F. 
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5.   TRANSFORMATION  OF  FIELD  EQUATIONS  FROM  A  FIXED  COORDINATE 
TO  A  MOVING  COORDINATE  SYSTEM 

As  discussed  previously,  the  reaction  rate  at  the  air  in- 
let surface  of  the  porous  medium  results  in  surface  recession. 
To  account  for  this  surface  recession,  the  field  equations 
must  be  transformed  from  a  fixed  coordinate  to  a  moving  coor- 
dinate system.   The  method  of  transforming  the  field  equations 
will  be  shown  for  only  the  porous  solid  heat  transfer  equa- 
tion since  the  approach  for  the  other  field  equations  (i.e., 
air  heat  transfer,  and  combined  Darcy ' s  law  and  continuity 
equation)  is  identical. 

The  porous  solid  heat  transfer  equation  with  the  modified 
heat  generation  term  is. 


ST 

|-[(l-p)  (k  +k  )-^]  -  hz(T  -T  )  +  R'6^(x=  0)  (A.  46) 

ax""    ^   e   r  9x         c   a     g 

8T 
=   ^l-P^^c^cTF 


Since  the  x  coordinate  is  a  function  of  time  during  surface 
recession  (e.g.,  T  (x(t),t)),  the  time  derivative  term  in 
equation  A. 4  6  must  be  expanded  using  the  chain  rule. 


3T     ,       3T 
'^     [T^(x(t),t)]     =   (^)  (§$)   +  (^) 


St'  c 


dx  '     'dt'  .        '   3t 

t     1        x 


(A. 47) 


The  X   coordinate  in  the  field  equations  is  nondimensionalized 
by  n  =  x/L.   Substituting  n  foi"  x,  expression  A. 47  becomes. 
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dx.dt.  ^^x  ^        ^^ti  X 


Noting   that. 


in     =     i 
9x  L 


(A. 49) 


and. 


dx      _      _    dn    _L        dL  dL  /tv    cn^ 

dt^^dt-'^dt^^dt  ^^-50) 


expression  A. 48  becomes. 


9T     ,        aT           ,_  9T     aT 
{ £)  {^)   +  ( £)    =   II(^) £  +  £  (A  51) 

^ax\^dt^   ^  at  ^     LMt^  an    at  {a.^^d 


Upon  substituting  n  =  x/L  into  the  remainder  of  equation  A. 46, 
and  substituting  expression  A. 51  for  the  time  derivative  term, 
the  porous  solid  heat  transfer  equation  for  the  surface 
recession  problem  is. 


aT 

l"^  i_[(l-p)  (k  +k  )_£]  -hz(T  -T  )  +R'6^(n=  0)  (A. 52) 

dr\  e   r  an        c   a    g 

.  aT   aT 

=   (l-p)p^c  [2-L^  +  3^] 

c  c  L   an    at 


As  Stated  in  Section  III.G,  the  thickness,  L,  as  a  function 
of  time,  and  L  can  be  obtained  from  expression  III.  44.  The 
air  temperature  and  combined  Darcy ' s  law  and  continuity 
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equation  appear  similar  to  equation  A. 52  for  the  surface 
recession  problem. 
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APPENDIX  B 
FORMULATION  OF  AUXILLIARY  EQUATIONS 

1.   CONTINUITY  EQUATION 

The  law  of  conservation  of  mass  requires  that  the  substan- 
tial derivative  of  the  fluid  mass  in  a  differential  volume  be 
zero.   Therefore,  the  continuity  equation  for  a  fluid  in  a 
porous  medium  is  expressed  by. 


^^(pp^dV)   =   0  (B.l) 


or  in  an  equivalent  form. 


I^(PP^)  +u|3^(pp^)  -PP,|^  =   0  (B.2) 


2.   RADIATION  HEAT  TRANSFER  ANALOGY  TO  FOURIER'S  LAW 

Radiation  heat  transfer  in  the  model  was  represented  by 
an  analog  to  Fourier's  law  of  conduction  heat  transfer  shown 
by. 


■Jrad  =  -"^r  ^  '^-2' 


where  k   is  an  equivalent  conductivity  of  the  particles  due 
to  radiation.   The  development  is  as  follows.   Assuming  air 
to  be  transparent  to  radiation,  and  treating  the  idealized 
geometry  of  the  porous  medium  (Figure  III.l)  as  a  series  of 
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closely  spaced  walls,  the  net  radiation  heat  flux  between 
adjacent  walls  is, 


[  .      =   ^-^(T^  -  T^^,  )  (B.4) 

^rad     2-e   x    x+dx  ' 


where  T  and  T   ,   are  the  respective  absolute  wall  tempera- 
tures, e  is  the  emissivity  of  the  carbon,  and  a  is  the  Stefan- 
Boltzman  constant.   Expanding  q   ,  in  a  Taylor  series  about 
T  .J,./  and  neglecting  higher  order  terms,  the  series  expan- 
sion may  be  written  as. 


'rad     2-£   x    x    2-e   x+dx  9x 


Simplifying,  the  above  expression  becomes. 


4ae  "3     9T  ,  ,^    ^. 

q   J   =   -  n T  ,  ,   TT—  dx  ( B .  6 ) 

^rad       2-e   x+dx  dx 


Equating  expressions  B.3  and  B.6,  and  noting  that  9T/3x  =  3T/3x, 

k  becomes, 
r 


'^r  =  5?f'^-W  <^-'' 


where  dx  is  now  equal  to  6,  the  pore  diameter.   From  the  close 
spacing  of  the  particles,  the  temperature  difference  will  be 
small  as  compared  to  the  magnitude  of  the  temperature.   Noting 
this,  the  average  absolute  temperature  of  the  carbon  particle 
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may  be  substituted  for  T   ,  .   The  equivalent  radiation  con- 
ductivity expression  becomes. 


k    =   ^P  (B.8) 

r     2-e    c 


and  the  radiation  heat  transfer  from  particle  to  particle 
may  be  represented  by. 


^rad        2-e   c  8x 


For  small  pore  diameters  equation  B.9  will  be  a  good  approxi- 
mation of  the  radiation  heat  transfer  between  particles. 

3.   POLYNOMIAL  APPROXIMATIONS  OF  THERMAL  PROPERTIES 

Relations  giving  the  dynamic  viscosity,  thermal  conduc- 
tivity, and  specific  heat  at  constant  pressure  of  air  at 
different  temperatures  were  required.   A  simple  method  to 
obtain  values  for  these  properties  is  to  fit  empirical  data 
with  2nd  order  Lagrange  polynomials. 

The  general  form  of  the  2nd  order  Lagrange  plynomial  is. 


(T.-T„)  (T.-T-)        (T.-T^)  (T.-T,  ) 
^1     (T^-T2)  (^^-^3)   1    (T2-T^)  (T2-T^)  ^2  ^n.JLV) 

(T.-T,) (T.-T„) 


(T.,-T,  )  (T.,-T^)   3 


where  k.  is  the  property  value  at  the  ith  temperature,  T.. 
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Choosing  three  temperatures  that  are  representative  of 
those  observed  during  the  analysis,  the  corresponding  values 
of  the  properties  are, 


sc 

:ript) 

Temp. 

y 

k 
a 

c 
a 

deg-F 

Ibm/ft-hr 

Btu/ft-hr-F 

Btu/lbm-F 

1 

80 

4.7973x10"^ 

1.5161x10"^ 

.24020 

2 

1340 

1.0045X10"-'- 

3.9013x10"^ 

.27268 

3 

3140 

1.5725x10"^ 

7.1647x10"^ 

.31196 

Applying  expression  B.IO  to  each  set  of  properties  results 
in  the  following  set  of  polynomials. 


y      =      -3.308x10"^      T^    +    4.633 xlo"^T     +4.427xl0~^  (B.ll) 

a  a 


k        =      -2.608  X  lO"-'-^   T^    +    1.930  xl0"^T      +    1.361xlo"^       (B.12) 
a  a  a 


c        =      -1.293x10"^      T^    +    2.758x10    ^T      +    .238  (B.13) 

a  a  a 


Each  expression  gives  property  values  within  two  percent  of 
the  data  for  temperatures  to  3000  degrees  Fahrenheit. 

4.   JUSTIFICATION  OF  THE  DANCKWERT ' S  BOUNDARY  CONDITIONS 
The  Danckwerts '  boundary  conditions  for  the  air  heat 
transfer  are. 


3T 
^a  TF  =   Pa  ^^a^^a-T«.)  ^^-l^^ 
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3T 


Bischoff  [41]  presents  a  discussion  of  these  boundary  condi- 
tions applied  to  mass  diffusion  equations  in  porous  media. 
An  analogous  discussion  is  made  here  for  the  Danckwerts' 
boundary  conditions  applied  to  the  air  heat  transfer  equation. 
The  porous  medium  with  entrance  and  exit  regions  is  shown 

in  Figure  B.l,   The  air  temperature,  T  ,    for  each  section  of 

a 

the  region  is  distinguished  by  subscripts,  as  are  the  proper- 
ties.  The  air  heat  transfer  equations  for  each  of  the  regions 
are  as  follows, 

d^T  dT 

a  a^ 

k   ^  -  p  u,  c   —5 =   0       x<0      (B.16) 

a   ,  2     ^a  1  a   dx  — 

dx 


,       dT  dT 

d^T  dT 

a_  a  .^ 

k   ^  -  p  u^  c   -g— ^  =0       X  >  L      (B.18) 

a   ^2  ^a  2  a   dx  - 


with  the  boundary  conditions. 


T   (-0°)   =   T  (B.19) 

^1 


T^  (0)   =   T^(0)  (B.20) 

a-.  a 
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Flow 

={> 

Entrance 
Region 

Po 

rous 

Slab 

Exit 
Region 

a^    a^  1 

V^ 

a'^ 

^2    ^2  2 

x  =  0 


x  =  L 


FIGURE  B.l  Justification  of  the  Dankwerts 

boundary  conditions. 
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Pa^l^a^^O^  -Ky\^'^^    =  PPa-^a^a^O^  -P^aS^tT^^O)] 


(B.21) 


T^  (L)   =   T^(L)  (B.22) 

^2  ^ 

Pa^2Va  (1)  -^l^t%(^)^  =  PPa^^a^a^Q)  -P^^[T  (L)  ] 
2  2  ax   a 

(B.23) 

T   (oo)   =   finite  (B.24) 

^2 


The  above  boundary  conditions  impose  the  restriction  that 
there  is  no  convection  heat  transfer  from  the  porous  solid 
to  the  air  at  the  boundary  surfaces.   The  difficulty  of  mixing 
Danckwerts'  and  convection  heat  transfer  boundary  conditions 
is  discussed  in  Section  III.H.   For  convenience,  the  thermo- 
physical  properties  will  be  treated  as  constant  in  the  entrance 
and  exit  regions. 

An  analytical  closed-form  solution  of  this  set  of  equations 
is  unknown  because  of  the  nonlinearity  of  the  temperature 
dependent  properties  in  equation  B.17.   However,  the  solutions 
of  equations  B.16  and  B.18  are. 


p  u,  c 
T^    =   K-|_  +  K2  exp  (  ^]^      X)  x  £  0      (B.25) 

1  a 


P  u^  c 
T^    =   K3  +  K^  exp  i    ^^  X)  X  ^  L      (B.26) 

2  a 
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Applying  boundary  conditions  B.19  and  B.24,  the  following 
results  are  obtained. 


K,   =   T  (B.27) 

1      °° 


K^   =   0  (B.28) 


The  solutions,  B.25  and  B.26,  become, 


p  u,  c 
T    =   T   +  K^  exp(  ^  , — -   x)  X  <  0      (B.29) 

1  a 


T    =   K.  X  >  L      (B.30) 

32       J 


It  would  be  necessary  to  have  the  solution  for  the  nonlinear 
equation  B.17  to  solve  for  K2  and  K^.   However,  the  constants 
need  not  be  known  to  continue  with  the  analysis. 
From  equation  B.29, 


T   (0)   =   T   +  K^  (B.31) 

a-,         «>    2 


and. 


|.[T   (0)]   =  '-^^^2  '2.32) 

1  a 


Substituting  these  into  equation  B.21  gives. 
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Pa^lVco^Pa^l^a^2-Pa^lV2  =  P^a^^a^a^^^  "  P^^i^^^a  ^^^  ^ 

(B.33) 

Cancelling  terms  and  noting  u,  =  up  (i.e.,  i^-iP-i  =  ^2^9^ 
yields, 


pp^uc^T^   =   PP^uc^T^(O)  -  pk^^[T^(0)]       (B.34) 


Rearranging,  the  Danckwerts '  boundary  condition  at  x  =  0  is. 


■^a  lrf'a<°'l   =  Pa"=a<^a  "  ^-'  '■^•^^'> 


Noting  that  K-.  is  a  constant,  the  derivative  of  equation  B.30 
is. 


^  T    =0  (B.36) 

dx  a^ 


and. 


|j[T^^(L)l   =   0  (B.37) 


Sbustituting  these  expressions  and  equation  B.22  into  equa- 
tion B.23,  and  noting  U2  =  up  (i.e.,  u,p-  =  U2P2)  gives. 


PPa'^'^a^a'"  -  ''a  aST'^^a'^*!   =   PPa'^'^a^a*"     '^•2^' 


Upon  cancelling  terms,  the  second  Danckwerts'  boundary 
condition  becomes, 
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^[T^(L)]   =   0  (B.39) 


An  important  consideration  for  using  the  Danckwerts ' 
boundary  conditions  is  that  it  simplifies  the  analysis  since 
equations  similar  to  B.17  may  be  solved  independently  without 
having  to  consider  entrance  and  exit  regions. 

5.   PRESSURE  DIFFERENTIAL  CALCULATION  RESULTING  FROM  EXTERNAL 
FLOW  AT  SURFACE  X  =  L 

To  obtain  the  pressure  differential,  AP,  for  simulating 

the  conditions  of  Fontenot's  experiments  [1]  (discussed  in 

Section  III.H),  Bernoulli's  equation  was  used.   The  following 

observations  showed  this  to  be  a  valid  assumption.   Schlichting 

[54]  points  out  that  for  the  ratio  of  u/U   in  the  range  of 

.0001  to  .01,  the  effects  of  "blowing"  or  "suction"  on  the 

potential  flow  over  the  external  surface  of  the  porous  medium 

may  be  neglected.   A  typical  value  of  u/U^  for  the  model  at 

which  U  =  25  knots  was  .0028.   For  steady  flow  over  a  flat 

plate,  the  flow  field  outside  the  boundary  layer  may  be 

described  by  Bernoulli's  equation.   This  is  a  direct  result 

of  the  Navier-Stokes  equation.   In  addition,  the  pressure 

gradient  across  the  boundary  layer  may  be  taken  as  zero. 

Therefore,  Bernoulli's  equation. 


2 
f  —  +  ^  =   constant  (B.40) 

Pa    2 


may  be  used  to  obtain  the  pressure  differential  across  the 
plate.   The  parameters  in  equation  B.40  are  defined  as  follows 
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P,  pressure;  p  ,  density  of  air;  U,  free  stream  velocity  of 
a 

the  air  over  an  external  surface  of  the  porous  medium.   For 
the  model,  the  density  of  the  air  was  approximated  by  the 
ideal  gas  law  as, 


p^   =   P/R^  T  (B.41) 

a.  a   a 


where  R   is  the  gas  constant  of  air,  T   is  the  absolute  tem- 
a        ^  a 

perature  of  the  air,  and  P  is  the  pressure.   Substituting 
equation  B.41  into  equation  B.40  gives. 


R   T       „2 
/   ^   ^  dP  +  ^  =   constant  (B.42) 


Upon  integrating,  equation  B.4  2  becomes. 


2 

U 
R   T   ln(P)  +  ^5-  =   constant  {B.43) 

a   a  ^ 


or, 

2  2 

U^         .  U^ 

R^  T    ln(P,)  +  -^  =   R   T    ln(P_)  +  -:f  (B.44) 

a   a->      J.     z       a   a.^      z     z 


Letting, 


P,  =  P  ,   U,  =  0,   T    =  T    =  T  ,   P„  =  P^  ,  U^  =  U 
1     »     1        a,     a^     «>'    2     L    2     = 

1      2 


and  substituting  these  into  equation  B.44,  yields. 
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R^    T   In(P^)   =   R^  T   ln(P_)  +  -^  (B.45) 


Rearranging,  equation  B.45  becomes, 

R      T   ln(P  /P^)   =  -^  (B.46) 


Taking  the  exponential  of  both  sides,  and  solving  for  P  , 

Li 


results  in. 


P^   =   P   exp(-U^/2R   T  )  (B.47) 

L       00    1^  ^   oo^    a   OO' 


From  the  above  expression,  and  noting  that  AP  =  P^  -  P^ ,  AP 


may  be  expressed  as. 


AP   =   P^[exp(-uf/2  R^  T  )  -  1]  (B.48) 

oo    *■    00    a  ^^ 


Expression  B.48  was  used  to  approximate  the  pressure  differ- 
ential across  the  porous  medium  when  simulating  the  conditions 
of  Fontenot's  experiments  [1]. 
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APPENDIX  C 
NUMERICAL  FORMULATION 

1.   FINITE  ELEMENT  METHOD 

The  solution  of  the  system  of  nonlinear,  coupled,  partial 
differential  equations  given  by  equations  III. 32,  III. 38,  and 
III.  39,  si±>ject  to  boundary  and  initial  conditions,  was  ob- 
tained by  a  Galerkin  formulation  of  the  finite  element  method, 
The  solution  of  equation  III. 12  was  obtained  using  a  shooting 
method . 

a.   Galerkin  Formulation 

A  Galerkin  formulation  of  the  Finite  Element  Method 
was  used  to  obtain  solutions  of  the  porous  solid  and  air 
energy  equations,  and  the  oxygen  diffusion  equation.   A  con- 
venient form  of  equations  III. 32,  III. 38,  and  III. 39  was 
used  in  the  formulation  as  shown  by, 

ST  9  T 

l"^   IrrLd-p)  (k   +k^)^]  -hz(T  -T^)  +R   z      =       (l-p)p    c^^ 
dx]  erdri  ca  g  ccdt 

(CD 


!^T  9  T 

l"2   |_(pk     _^)   -pp^  c     uL"1    ^+hz(T^-T^)  (C.2) 

dri  ^    on  a     a  dri  c      a 

3  T 
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where  the  spatial  coordinate,  x,  was  nondimensionalized  by 
n  =  x/L. 

The  closed  domain  (0,1)  was  partitioned  into  (n-1) 
contiguous  elements  of  variable  length  i.,    i  =  l,...,n-l. 
This  defines  an  n  nodal  point  model.   The  three  field  varia- 
bles, T    ,    T    ,    <p   were  approximated  by, 
c   a 

T^(n,t)  =  ij;^(n,t)  =  I   G(n)9^(t)  (c.4) 

T^(n,t)  =  ^^{n,t)     =  [  G(n)e2(t)  (c.5) 

<l>(n/t)  =  H)^{T],t)     =     I   G(n)e^(t)  (c.6) 


where  G.,  for  i  =  l,...,n  is  a  set  of  specified  basis  func- 
tions with  local  support,  and  the  sets  {Q^,Qy,Q~,;    i  =  l,...,n} 

are  the  solution  coefficients  to  be  determined.   The  G.  were 

1 

]^ 
selected  to  atisf y  the  condition  G. (n-)  =  5. .,  where  the 

k  k 

Kronecker  delta,  5. .,  is  defined  by  5 . .  =  1  for  i  =  j,  and 

]^ 

5..  =  0  for  i  7^  j  .   As  a  result,  9,,  9-,  and  9-,  are  the 

values  i|^,  ,  i)y,    ij;,  at  the  nodal  points  (i.e.,  9,  (t)  =  i|;,(ri-/t)) 

i 
Linear  interpolation  functions  (shown  in  Figure  C.l) 

were  used  as  the  basis  functions.   These  are  the  lowest 

polynomial  functions  which  provide  the  necessary  function 

continuity. 

As  a  measure  of  error,  a  residual  function,  r.,  is 

defined  for  each  field  equation  by, 
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n-2   n-1 


n 


FIGURE  C.l  Linear  shape  functions  used  in  the 

Galerkin  formulation  of  the  FEM. 
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r^(n,t)   =  K^ii,)    -  ^^        i  =  1,2,3  (C.7) 

where  A.  denotes  the  spatial  operator  of  the  ith  equation. 
For  convenience,  the  following  convention  for  differentiation 
is  adopted, 


!-()=()'  (C.8) 


.2 

^(   )   =   (   )"  (C.9) 

3n 


|^(   )   =   (  *  )  (C.IO) 


For  field  equations  C.l,  C.2,  and  C.3,  the  residuals  are. 


-  n  n 

r,   =   l"^[(1-p)  (k  +k  )  J  G!9.  ]'-h2  I      G.  (9^  -9^  )  (C.ll) 
-L  ®   ^  i=l  ^  ^1       i=l   ^   -^i   ^i 


,    n  n     . 

+  Re"  z  I    G  9    -  (l-p)p  c    I      G  9 
^  "^i   i=l  ^  -^i         ^  ^  i=l   ^  -"i 


-9      '^  -1   ^ 

r^   =   L   (pk^  y  G!9o  )'  -  pp^c^uL  ^      I      G!9„       (C.12) 

^  ^  i=l  ^  "^i        ^  ^       i=l   ^  ^i 


n  n     . 

+  hz   j;   G.  (9,  -9,  )  +  uL"^(pP)  '  -  p  p^  c    I      G.9^ 
i=l   ^   -^i   ^i  ^  ^  i=l   "■  ^i 
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r3   -   L"^(P^e  I   G|e   )•  -  l"^  (p  u  f     G.9   )•        (C.13) 

i=l     i  i=l     i 


_,    n  n 


-  «o  83:^.^,  V3.  -  P  J,  Vs. 

2   1   1=1      1      1=1      1 


where  the  coefficients  multiplying  the  response  variables  are 
themselves  functions  of  the  response  variables,  and  thus, 
the  equations  are  nonlinear.   In  accordance  with  the  Galerkin 
method,  the  final  system  of  ordinary  differential  equations 
was  obtained  by  setting  each  residual,  r.,  orthogonal  to  each 
basis  function,  G.,  that  is. 


1  i=l,2,...,n 

/   G.  r .  dn   =   0  (C.14) 

I       -■  j  =  1,2,3 


The  3n  ordinary  differential  equations  given  by  equations  C.14 
retain  the  character  of  the  original  set  of  partial  differ- 
ential equations.   Thus,  linear  field  operators  transform  to 
matrix  operators  and  nonlinear,  coupled  field  operators  become 
nonlinear,  coupled  algebraic  operators.   Incorporation  of 
the  boundary  conditions  resulted  in  3n  nonlinear  coupled 
ordinary  differential  equations, 

A(t)  ^{B^,Q^,Q^)    +   F(t)   =   BCt)  ^p  (C.15) 

subject  to  initial  conditions,  where  B  is  a  3n x  3n  matrix. 
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A  is  the  operator  associated  with  the  field  operator  A .  in 
expression  C.l ,    and  F  is  an  excitation  vector. 
Adopting  the  convention. 


n 
<G.>  {9.}   =    I      G.9.  (C.16) 

1     1       i=l   ^  ^ 


and  applying  the  operation  of  expression  C.14  with  an  integra- 
tion by  parts  on  the  second  order  derivatives  gives. 


{G^}l"^(1-P) (k^+k^)<G^>'{e^}|J   -   l"^(1-p) (k^+k^)    /    {G^} •<G^> 'dn 


1  I 

-hz      /    {G. }<G.>dn{9,}    +   hz      /    {G. }<G . >dn (9    }  (C.17) 

0^3  J-  0^^ 


,1  1  , 

+    R   z9,         /    {G.  }<G.>dn{9^}    =    (l-p)p    c         /    {G  .  }  <G  .  >dn  (  9  ,  } 
gjQ^i:)  3  ^^^cc    Q^ij  1 


{G.}L~^pk^<G,>'{9-,}|i    -    L    ^pk^      /    {G  .  }  '  <G  .  > 'dn  {9„  }  (C.18) 

1  a      J.  z      vj  a    ^  13  z. 


_-,         1 

pp    c    uL  /    {G.  }<G.>'dn{9»}    +   hz  /{G.  }<G.>dn{9,  } 

aa^i;]  z  1;)  1 


1  -1  ^ 

■hz      /    {G.  }<G.>dn{9^}    +   uL    "^(pP)'       /    {G.}dn 

0  ^         ^  0  ^ 


1 

P  P^  C;,      /    ^^i  }<G    >dn{9^} 
a     a   Q^         13  / 
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{G.}L"^pf    <G.> '{9    }  1^    -   L~^pl?         /    (G.} '<G.>'dn{9    }  (C.19) 


-1      1  1  1 

puL    ""      /    {G.}<G.>'dn{9    }    -    L       (pu)  •      /    {G .  }<G .  >dn{9 .} 
0^  -^  O'"'^ 


-1     1  1  • 

R^  29^    /  {G.}<G.>dn{9  }   =   p   /  {G.}<G.>dn{9  } 

^2    0^^    J      qi:i    j 


The  first  term  in  each  of  the  above  expressions  is  a  boundary 
term  which  permits  incorporation  of  natural  boundary  condi- 
tions.  Implementation  of  the  boundary  conditions  is  presented 
in  Section  C.l.b.   The  coefficients  in  equations  C.17,  C.18, 
and  C.19  are  comprised  of  variable  dependent  properties,  and 
were  taken  as  the  average  value  of  the  properties  over  an 
element.   In  the  limit,  as  the  elements  get  smaller  (i.e., 
n  ->■  <»)  ,  the  average  values  of  the  coefficients  converge  to 
the  exact  values. 

Inspection  of  expressions  C.17,  C.18  and  C.19  shows 
the  four  operators, 

/  {G.}'<G.>'dn  (C.20) 

/  {G.}<G. >'dn  (C.21) 

/  {G.}<G. >dn  (C.22) 

/  {G^}dn  (C.23) 
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To  formulate  these  operators,  the  global  shape  function, 
G.,    was  defined  on  the  local  level  by. 


(i-1)  ^  ^(i) 


i  =  g|^-^^©g^^^  (C.24) 


where  g,  and  g^  were  defined  by. 


(1  -  Y~)  for  E,    in  element  (e) 


e 


g{®^   =  (C.25) 

0        for  E,    not  in  element  (e) 


for  E,    in  element  (e) 


g2®^   =  (C.26) 

0  for  E,    not  in  element  (e) 


and  I      is  the  length  of  the  eth  element.   The  (+)  notati 


e 


on 


in  expression  C.24  means  that  G.  is  the  union  of  g,      and 
g^"^  .   The  local  shape  functions  (i.e.,  the  elements)  have 
the  following  properties, 


Z 
e 

(i)       /    g|^^  gi"^^      =   0     if   j  ?^  m  (C.27) 

0     IK 


1   if  i  =  j 
(ii)     g/^^  (n J  =  6^.  =  (C.28) 


1    2  i: 


0   if  i  7^  j 
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Having  defined  the  local  shape  functions,  the  ele- 
mental matrix  operators  contributing  to  the  global  matrix 
operators  C.20  through  C.23  are. 


/  {G.}'<G  >'dn  -  ^ 


-1 


-r 


(C.29) 


/  {G^}<G.>' 


-1 


Ij 


(C.30) 


/  ^G.}<G  >   -   -^ 


_1 


(C.31) 


/  ^G.} 


(C.32) 


The  derivations  of  these  operators  are  presented  in  Section 
C.l.d  of  this  appendix. 

b.   Implementation  of  Boundary  Conditions 

Having  formulated  the  system  matrices  for  the  field 
equations,  treatment  of  the  boundary  conditions  will  now  be 
discussed.   Each  field  equation  is  considered  individually. 
1.   Porous  Solid  Transfer  Equation 

The  third  set  of  porous  solid  heat  transfer 
boundary  conditions  (i.e.,  expressions  III. 57  and  III. 58)  will 
only  be  considered  since  the  first  and  second  sets  are  sub- 
sets of  the  third.   The  third  set  of  boundary  conditions  for 
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the  porous  solid  are. 


9  Q 

L"-'-(k  +k  )-^  =   h^  (T  -T  )  +  ae{T^-T^)         n  =  0        (C.33) 
e   r  8ri       1   c  o°       '  c   <» 

38 
L"-^(k  +k  )^  =   -h_(T  -T  )  -  ozii^-il)       n  =  1        (C.34) 


Since  the  first  term  in  expression  C.17  is. 


-2  ^ 

{G.}L  ^(l-p)(k  +k  )<G.>'{9,}|  (C.35) 

1  e   r    1     1   Q 


or  in  analogous  form, 

3  6 
l"^(1-P) (k^+k^)  -j^  (C.36) 

Natural  boundary  conditions,  C.3  3  and  C.34,  may  be  directly 
substituted  into  equation  C.36.   The  response  dependent 
parameters,  h, ,  h« ,  and  T  ,  changing  with  time,  are  evaluated 
continuously  within  the  integration  routine.   Thus,  the  boun- 
dary conditions  are  incorporated  in  the  system  matrices  as 
follows . 


(1)  -L   (l-p)h, :   added  to  the  stiffness  matrix  A(t) 

at  location  A,  ,  ~ 

(2)  l"-'- (1-p)  [h^T  -a£  (t'^-t'*)  ]  :   added  to  the  excitation 

vector  F(t)  at  location  F, 
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(3)  -L   (l-p)h^:   added  to  the  stiffness  matrix  A(t) 

at  location  A-,   ^  -,   ^        ~ 
3n-2 , 3n-2 

(4)  l"-'-(1-p)  [h^T  -ae  (T  -t"^)]:   added  to  the  excitation 

2°°         Coo 

vector  F(t)  at  location 
P 


2.   Air  Heat  Transfer  Equation 

The  first  and  second  set  of  boundary  conditions 
for  the  air  heat  transfer  equation  are. 


3T 

l"-"- k  ^^  =   p  c  u(T  -T  )  n  =  0      (C.37) 

a  3n      ^a  a    a  00 

3T 

-5-^=0  n  =  1      (C.38) 

3n 


Since  these  are  both  natural  boundary  conditions,  they  are 
substututed  for  the  first  term  of  expression  C.18  at  n  =  0 
and  n  =  1/  respectively.   The  time  dependent  properties  and 
parameters  in  the  coefficient  are  evaluated  continuously 
within  the  integration  routine.   The  boundary  conditions  are 
implemented  by  adding. 


(1)  -p p  c  u L   :   to  the  stiffness  matrix  A(t)  at  location 

(2)  pp  c  uL   T  :   to  the  excitation  vector  F(t)  at 

^       "*    location  F- 


The  third  set  of  air  heat  transfer  boundary  conditions  are. 


T   =   T  n  =  0      (C.39) 

a      <» 
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8T 

^  =  0  n  =  1     (C.40) 


The  essential  boundary  condition  at  n  =  0  is  imposed  in  the 
Galerkin  equation  as  follows.   The  A,,  .  row  of  the  A(t) 
matrix,  the  B^  .  row  and  the  B.  ^  column  of  the  B(t)  matrix, 
and  the  F-  location  of  the  excitation  vector,  F(t),  are  all 
set  equal  to  zero.   The  B2  2  location  of  the  B(t)  matrix 
is  then  set  equal  to  one. 

3.   Oxygen  Transport  Equation 

For  the  oxygen  diffusion  equation,  the  first  and 
second  set  of  boundary  conditions  are. 


l""^  ^o  I^  =  ^(<J^  -  't>J  n  =  0     (C.41) 


1^  =   0  n  =  1      (C.42) 

dn 


since  these  are  natural  boundary  conditions,  they  were  sub- 
stituted for  the  first  term  in  expression  C.19.   The  proper- 
ties are  evaluated  continuously  within  the  integration  routine 
The  boundary  conditions  are  implemented  by  adding. 


(1)  -p  u L   :   to  the  stiffness  matrix  A(t)  at  location 

A  ~ 

^3,3 

(2)  puL   (|)  :   to  the  excitation  vector  F(t)  at  location 

°     F3 
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The  third  set  of  boundary  conditions  for  the  oxygen  diffusion 
equation  are, 

^   =   <1)^  n  =  0      (C.43) 


1^  =   0  n  =  1      (C.44) 


The  essential  boundary  condition  at  n  =  0  is  imposed  in  the 
Galerkin  equation  as  follows.   The  A^  .  row  of  the  A(t) 
matrix,  the  B-.  .  row  and  the  B.  ^  column  of  the  B(t)  matrix, 
and  the  F^  location  of  the  excitation  vector,  F(t),  are  all 
set  equal  to  zero.   The  B^  ^  location  of  the  B(t)  matrix  is 
then  set  equal  to  one. 

4.   Surface  Recession  Problem 

As  discussed  previously,  the  oxygen  diffusion 
equation  is  eliminated  during  the  surface  recession  phase  of 
a  problem.   By  doing  this,  a  reordering  of  the  system  matrices 
must  take  place,  including  the  locations  for  applying  the 
boundary  conditions.   The  boundary  conditions  of  the  particle 
and  air  temperature  equations  are  implemented  the  same  as  above 
except  the  indices  identifying  the  matrix  locations  of  the 
form,  3n-i,  are  changed  to  2n-(i-l). 

c.   Treatment  of  the  Reaction  Rate  Term 

An  exponential  reaction  rate  term  appears  in  both  the 
porous  solid  heat  transfer  and  the  oxygen  diffusion  equation. 
The  modified  Gear  integration  routine  requires  calculating  or 
approximating  the  Jacobian  matrix  from  the  system  matrices , 
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A(t)  and  BCt),   In  order  to  include  the  reaction  rate  term 
in  the  Jacobian  matrix,  thus  improving  the  efficiency  of  the 
integration  routine,  the  exponential  terms  were  treated  as 
follows.   The  reaction  rate  terms  were  multiplied  and  divided 
by  the  oxygen  concentration  (only  the  heat  generation  term 
is  shown,  oxygen  consumption  is  treated  identically) , 


-1   ^ 
R  z9->    y   G.e.  (C.45) 

g    3 .  > ,   1  3  • 

^     1  1=1      1 


Applying  the  operation  of  expression  C.14,  and  moving  the 
reaction  rate  term  divided  by  the  oxygen  concentration  out- 
side the  integral,  the  Galerkin  operator  for  the  heat 
generation  becomes. 


R  z  e"-"-   /  {G.  }<G.>dn{6-,}  CC.46) 

g   3  Q^   1   :    3 


The  coefficient  of  expression  C.46  is  treated  as  a  time 
dependent  parameter.   The  Galerkin  operators  are  distributed 
into  the  stiffness  matrix  at  locations,  3n-2,3n  for  heat 
generation,  and  3n,3n  for  oxygen  consumption.   During  the 
surface  recession  phase,  reaction  rate  is  treated  in  a 
different  manner  as  discussed  in  Section  III.G. 
d.   Derivation  of  the  FEM  Operators 

In  the  section  on  the  finite  element  formulation, 
the  following  four  differential  operators  were  identified,  - 
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1 

/  {G?}<G!>dn  (C.47) 

0     ^    ^ 


1 
/  {G.}<G:>dn  (C.48) 

0     ^    ^ 


1 

/  {G.}<G.>dn  (C.49) 


1 

/  {G^}dn  (C.50) 


where  the  G.  are  the  global  basis  functions.   These  operators 
are  constructed  on  the  element  level  by  introducing  the 
corresponding  element  basis  functions,  g..   The  global  and 
element  basis  functions  are  related  by, 

G.   =  g{^-^'  ®  g^^*  (C.51) 

where  g,  and  g„  are  defined  by 


(1  -  •^)  for  ?  in  element  (e) 

g£^^   =  (C.52) 

0        for  E,   not  in  element  (e) 
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^2 


for  ^  in  element  (e) 


for  C  not  in  element  (e) 


(C.53) 


and  Z      is  the  length  of  the  (e)th  element, 
e  ^ 

The  derivation  of  the  local  elemental  matrices 
according  to  the  Galerkin  method  for  the  global  operations 
proceeds  as  follows: 
For  operator  C.47, 


global 


local 


/    {G'}<G!>dii 


0 


i^\-l 


<g^   g^>d^]       (C.54) 


Noting   that , 


(C.55) 


^2      =      — 


(C.56) 


the  elemental  matrix  becomes 


e 


-(i)2] 


i 


e 


dC 


-1 


-I 


(C.57) 
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For  operator  C.48, 


global 


local 


I 


1   {G.}<G:>dn 
0    ^   ^ 


0 


<g-^  g^>dC]  (C.58) 


Substituting  in  the  local  shape  functions  gives. 


(1  -  ^) 


e 


e   e 


(C.59) 


Carrying  out  the  operations  gives 


e        e 


.,4. 


2Z' 


1  .  ^  ^- 

^e  Z  ^ 
e 

2£'^ 


d^ 


-1    1 


_-l    1_ 


(C.60) 


For  operator  C.49, 


global 


local 


/  {G.}<G.>dn 
0     ^    ^ 


£ 


0 


g-^  g2>dC]  (C.61) 
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Substituting  in  for  the  local  shape  functions 


(1  - 


e 


<(1  - 


f)(f)>d^ 
e   e 


(C.62) 


the  elemental  matrix  becomes 


(1-2  — +  — 

e  i^ 


4    .2' 

e  i 

e 


4    ,2' 
e   £ 
e 


(1-2  f +M 

e   £ 
e 


dC 


1    =r 


L^   1 


(C.63) 


For  operator  C.50, 


global 


local 


/  (G^ldn 


0 


'e  Wl, 


d^] 


(C.64) 


Substituting  in  the  local  shape  function,  and  integrating, 
the  expression  becomes 


^e 


dC 


(C.65) 
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This  last  operator  is  used  for  the  excitation  vector  as 
described  in  the  FEM  formulation. 

2.   SHOOTING  METHOD 

The  solution  of  equation  III. 12  is  obtained  by  the 
shooting  method,  a  general  discussion  of  the  shooting  method 
with  examples  is  presented  by  Gerald  [55] .   The  method  is 
based  on  iterative  solutions  of  equation  III. 12  until  a  solu- 
tion is  reached  which  also  satisfies  the  boundary  conditions, 
P(0)  =  P^  and  P(l)  =  P^.   Along  with  the  boundary  condition 
at  n  =0,  an  initial  estimate  to  dP/dn  is  specified,  and 
equation  III. 12  is  integrated  from  n  =  0  to  n  =  1  using 
Euler's  method.   An  approximate  Newton's  method  is  then  used 
as  follows. 


dp     dp 
dP  dn  .  "  dn 

^\.  =  ill,.,  -^^^^i-iP(i),  -p(i);_,    ^^-^'^ 


dP 


i-2 


to  provide  a  better  approximation  of  dP/dn .   The  procedure  is 
repeated  until  the  solution  has  converged.   A  solution  satis- 
fying equation  III. 12  and  its  associated  boundary  conditions 
is  calculated  at  each  time  step.   The  current  values  of  the 

properties  are  used  and  9p  /3t  is  approximated  by  linear 

a 

I' 

interpolation  using  the  current  and  past  values  of  p  .   The 

a 

dp   /dt   term  was  neglected  for  the  first  two  integrations  of  a 
a 

problem.   The  reason  for  this  is  associated  with  the  difficulty 
of  specifying  a  reasonable  set  of  initial  conditions  for  the 
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temperatiore  and  oxygen  diffusion  equations.   When  equation 
III. 12  was  treated  as  an  initial  condition  problem,  the 
pressure  gradient  in  a  small  region  of  the  porous  medium  would 
approach  zero  for  some  initial  conditions  of  temperature  and 
oxygen  concentration  (the  starting  heat  generation  rate  for 
those  initial  conditions  was  large  at  time,  t  =  0) .   The 
pressure  gradient  approaching  zero  resulted  in  numerical 
instabilities  for  the  temperature  and  oxygen  diffusion  equations 
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